folds={'dycomsNewSubs','isdacNew2'};
ccn_clean={'200','100'};
divls=[3.75e-6,0];
titles={'DYCOMS','ISDAC'};
cdiams=[20 12; 18 10];

%Linewidths, symbols, and colors
lwinv=[3 3 3 1.5 .5 .5 3 .5];
symsinv='x_|pv<+o';
lwinvhr=[3 3 3 1.5 .5];
symsinvhr='x_|po';
cols=[1 .5 .5; .5 .5 1; .8 0 0 ; 0 0 .8];
cols2=[.9 .25 .25; .25 .25 .9];

%Legend text
lgtextinv={'None (M000)','Evap. Only (M100)','Sed. Only (M010)',...
    'Rad. Only (M001)','Sed. + Rad. (M101)','Evap. + Rad. (M011)','Sed. + Evap. (M110)','All (M111)'};
lgtitle={'{\color[rgb]{1,.5,.5}Clean}/{\color[rgb]{.5,.5,1}Clean}',...
    '{\color[rgb]{.8,0,0}Polluted}/{\color[rgb]{0,0,.8}Polluted}'};
lgtitle2={'{\color[rgb]{.9,.25,.25}DYCOMS}/{\color[rgb]{.25,.25,.9}ISDAC}'};
lgtextpolarinv={'None','Sed. Only','Sed.+Rad.','Rad. Only','Evap.+Rad.','Evap. Only','Sed.+Evap.','Sed. Only','All'};
lgtextdiff={['Evap. Only',10,'(C\_M100P - C\_M000)'],['Sed. Only',10,'(C\_M010P - C\_M000)'],...
    ['Rad. Only',10,'(C\_M010P - C\_M000)'],['All',10,'(C\_M100P - C\_M000)'],...
    ['Evap. Only',10,'(P\_M000 - P\_M100C)'],['Sed. Only',10,'(P\_M000 - P\_M010C)'],...
    ['Rad. Only',10,'(P\_M000 - P\_M010C)'],['All',10,'(P\_M000 - P\_M100C)']};
alphabet='a':'z';

%panel positioning
apos=[.1159 .5662 .2846 .3588; ...
    .4701 .5662 .2846 .3588; ...
    .1159 .11 .2846 .3588;...
    .4701 .11 .2846 .3588];
%% Make smaller file with averaged data

%Can only run this section if there is access to the full simulation output

dz=10;%meters
nz2=300;
CTind=150;
CTz=([1:300]-CTind)*dz;
ctd=-60:10:-10;
[we,radcool,mass,mLWP,entheat]=deal(nan(8,2,length(folds),2));
[radtop,evaptop]=deal(nan(8,2,length(folds),2,length(ctd)));
[mFTH,mRCP,mVAP,mRTP,mCCP,mTH]=deal(nan(nz2,8,2,length(folds),2));

for cd=1:2
    for ff=1:length(folds)
        if ff==1
            times=12:3:72;
        else
            times=12:3:108;
        end
        ntimes=length(times);

        %five number code interpretation:
        %ised ievp irad isat icldsed
        folders=cell(8,2,length(folds));

        cdiam=cdiams(ff,cd);

        folders{1,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_00001_0']; %full physics
        folders{1,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_00001_0'];
        folders{2,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_02001_',num2str(cdiam)]; %noevap
        folders{2,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_02001_',num2str(cdiam)];
        folders{3,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_10001_',num2str(cdiam)]; %nosed
        folders{3,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_10001_',num2str(cdiam)];
        folders{4,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_00101_',num2str(cdiam)]; %norad
        folders{4,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_00101_',num2str(cdiam)];
        folders{5,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_10101_',num2str(cdiam)]; %evap only
        folders{5,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_10101_',num2str(cdiam)];
        folders{6,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_02101_',num2str(cdiam)]; %sed only
        folders{6,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_02101_',num2str(cdiam)];
        folders{7,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_12001_',num2str(cdiam)]; %rad only
        folders{7,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_12001_',num2str(cdiam)];
        folders{8,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_12101_',num2str(cdiam)]; %deny all
        folders{8,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_12101_',num2str(cdiam)];

        for i=1:8
            for j=1:2
                display([cd,ff,i,j])
                %entrainment rate calculations - doing two ways
                rams_hdf5c({'theta','rv','rcp','pi','dn0'},[times(1),times(end)],folders{i,j,ff});

                runs.tempk=temp(runs.theta,runs.pi);
                runs.theta=runs.theta.*(1+0.61*runs.rv-runs.rcp);

                th1=squeeze(mean(mean(runs.theta(:,:,:,1))));

                newz=1:.01:length(th1);

                th1h=csaps(1:length(th1),th1,0.5,newz);
                th1hg=gradient(th1h);
                [~,indz]=max(th1hg);

                th1hg=th1hg/.01*100;
                indt=find(th1hg>0.5,1,'first');
                thbase=th1h(indt);
                indz2=find(th1h>thbase+2,1,'first');

                z1i=newz(indz)*dz; %* grid spacing
                z1i2=newz(indz2)*dz; %* grid spacing

                th1=squeeze(mean(mean(runs.theta(:,:,:,end))));

                th1h=csaps(1:length(th1),th1,0.5,newz);
                th1hg=gradient(th1h);
                [~,indz]=max(th1hg);

                th1hg=th1hg/.01*100;
                indt=find(th1hg>0.5,1,'first');
                thbase=th1h(indt);
                indz2=find(th1h>thbase+2,1,'first');

                z2i=newz(indz)*dz;
                z2i2=newz(indz2)*dz;

                dt = (times(end)-times(1))*10*60;
                we(i,j,ff,cd)=(z2i-z1i)/dt + (z2i+z1i)/2*divls(ff); %BL height change - average subsidence rate
                we2(i,j,ff,cd)=(z2i2-z1i2)/dt + (z2i2+z1i2)/2*divls(ff);

                %now calculate entrainment heating
                bltind=round(indz/100);
                meanT=squeeze(mean(mean(runs.tempk(:,:,:,1))));
                blT=trapz(meanT(1:bltind).*squeeze(runs.dn0(1,1,1:bltind,1)))/trapz(runs.dn0(1,1,1:bltind,1));
                dT=max(meanT)-blT;
                dntop=runs.dn0(1,1,bltind);

                entheat(i,j,ff,cd)=1004*we(i,j,ff,cd)*dT*dntop;

                %BL-top normalized profiles and integrals
                nz=size(runs.theta,3);

                nFTH = nan(nz2,ntimes);
                nVAP = nan(nz2,ntimes);
                nRCP = nan(nz2,ntimes);
                nRTP = nan(nz2,ntimes);
                nCCP = nan(nz2,ntimes);
                nDN = nan(nz2,ntimes);
                nTH = nan(nz2,ntimes);

                rams_hdf5c({'fthrd','vapliqt','nuccldrt','rcp','rrp','dn0',...
                    'theta','wp','ccp'},times,folders{i,j,ff});
                for t=1:ntimes
                    %find BL top
                    th1h=squeeze(mean(mean(runs.theta(:,:,:,t))));
                    th1h=gradient(th1h);
                    [~,indz]=max(th1h);

                    bottom = CTind - indz + 1;
                    top = bottom - 1 + nz;
                    nFTH(bottom:top,t) = squeeze(mean(mean(runs.fthrd(:,:,:,t))));
                    nVAP(bottom:top,t) = squeeze(mean(mean(runs.vapliqt(:,:,:,t)+...
                        runs.nuccldrt(:,:,:,t))))/600;%convert to kg/kg/s
                    nRCP(bottom:top,t) = squeeze(mean(mean(runs.rcp(:,:,:,t))));
                    nCCP(bottom:top,t) = squeeze(mean(mean(runs.ccp(:,:,:,t))));
                    nRTP(bottom:top,t) = squeeze(mean(mean(runs.rcp(:,:,:,t)+...
                        runs.rrp(:,:,:,t))));
                    nDN(bottom:top,t) = squeeze(mean(mean(runs.dn0(:,:,:,t))));
                    nTH(bottom:top,t) = squeeze(mean(mean(runs.theta(:,:,:,t))));
                end

                mnDN=mean(nDN,2);
                mnFTH=mean(nFTH.*nDN,2);
                radcool(i,j,ff,cd)=nansum(mnFTH(CTz<=0))*1004*dz;
                mnVAP=mean(nVAP.*nDN,2);

                for ic=1:length(ctd)
                    evaptop(i,j,ff,cd,ic)=sum(mnVAP(CTz>=ctd(ic)&CTz<=10))*2.5e6*dz;
                    radtop(i,j,ff,cd,ic)=sum(mnFTH(CTz>=ctd(ic)&CTz<=10))*1004*dz;
                end

                mLWP(i,j,ff,cd)=nansum(mean(nRTP.*nDN,2))*dz;
                mFTH(:,i,j,ff,cd)=mean(nFTH,2)*3600; %K/hr
                mRCP(:,i,j,ff,cd)=mean(nRCP,2);
                mCCP(:,i,j,ff,cd)=mean(nCCP,2);
                mRTP(:,i,j,ff,cd)=mean(nRTP,2);
                mVAP(:,i,j,ff,cd)=mean(nVAP,2)*2.5e6/1004*3600; %K/hr
                mTH(:,i,j,ff,cd)=mean(nTH,2);

            end
        end
    end
end
save('isdacNew_dycomsNewSubs5.mat','we','we2','evaptop','radtop','radcool',...
    'mLWP','mFTH','mRCP','mVAP','mCCP','mRTP','mTH','entheat')

%% Now do average files for high resolution simulations
folds={'dycomsNewSubs_hr','isdacNew2_hr'};

%using a different array structure than for the low resolution
%simulations... very annoying in retrospect

dz=3;%meters, stretched above 1km, but checked and can confirm that the BL ...
%top never exceeds 1km, so do not need to account for the stretching
nz2=900;
CTind=450;
CTz=([1:900]-CTind)*dz;
ctd=[-60,-51,-39,-30,-21,-9];
[we,we2,radcool,mass,mLWP,entheat]=deal(nan(10,length(folds)));
[radtop,evaptop]=deal(nan(10,length(folds),length(ctd)));
[mFTH,mRCP,mVAP,mRTP,mCCP,mTH]=deal(nan(nz2,10,length(folds)));
folders=cell(10,length(folds));
figure;hold on
for ff=1:length(folds)
    %output for high res simulations is every 15 minutes. It was every
    %10 minutes for the low res simulations
    if ff==1
        times=8:2:48;
    else
        times=8:4:72;
    end
    ntimes=length(times);

    %five number code interpretation:
    %ised ievp irad isat icldsed

    folders{1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_00001_0']; %full physics
    folders{2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_02001_',num2str(cdiams(ff,2))]; %noevap
    folders{3,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_10001_',num2str(cdiams(ff,2))]; %nosed
    folders{4,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_00101_',num2str(cdiams(ff,2))]; %norad
    folders{5,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_12101_',num2str(cdiams(ff,2))]; %deny all
    folders{6,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_00001_0'];
    folders{7,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_02001_',num2str(cdiams(ff,1))]; %noevap
    folders{8,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_10001_',num2str(cdiams(ff,1))]; %nosed
    folders{9,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_00101_',num2str(cdiams(ff,1))]; %norad
    folders{10,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_12101_',num2str(cdiams(ff,1))]; %deny all

    for i=1:10
        [ff,i]
        rams_hdf5c({'theta','rv','rcp','pi','dn0'},[times(1),times(end)],folders{i,ff});

        runs.tempk=temp(runs.theta,runs.pi);
        runs.theta=runs.theta.*(1+0.61*runs.rv-runs.rcp);

        th1=squeeze(mean(mean(runs.theta(:,:,:,1))));

        newz=1:.01:length(th1);

        th1h=csaps(1:length(th1),th1,0.5,newz);
        th1hg=gradient(th1h);
        [~,indz]=max(th1hg);

        th1hg=th1hg/.01*100;
        indt=find(th1hg>0.5,1,'first');
        thbase=th1h(indt);
        indz2=find(th1h>thbase+2,1,'first');

        z1i=newz(indz)*dz; %* grid spacing
        z1i2=newz(indz2)*dz; %* grid spacing

        th1=squeeze(mean(mean(runs.theta(:,:,:,end))));

        th1h=csaps(1:length(th1),th1,0.5,newz);
        th1hg=gradient(th1h);
        [~,indz]=max(th1hg);

        th1hg=th1hg/.01*100;
        indt=find(th1hg>0.5,1,'first');
        thbase=th1h(indt);
        indz2=find(th1h>thbase+2,1,'first');

        z2i=newz(indz)*dz;
        z2i2=newz(indz2)*dz;

        dt = (times(end)-times(1))*15*60;
        we(i,ff)=(z2i-z1i)/dt + (z2i+z1i)/2*divls(ff); %BL height change - average subsidence rate
        we2(i,ff)=(z2i2-z1i2)/dt + (z2i2+z1i2)/2*divls(ff);

        bltind=round(indz/100);
        meanT=squeeze(mean(mean(runs.tempk(:,:,:,1))));
        blT=trapz(meanT(1:bltind).*squeeze(runs.dn0(1,1,1:bltind,1)))/trapz(runs.dn0(1,1,1:bltind,1));
        dT=max(meanT)-blT;
        dntop=runs.dn0(1,1,bltind);

        entheat(i,ff)=1004*we(i,ff)*dT*dntop;

        %BL-top normalized profiles and integrals
        nz=size(runs.theta,3);

        nFTH = nan(nz2,ntimes);
        nVAP = nan(nz2,ntimes);
        nRCP = nan(nz2,ntimes);
        nRTP = nan(nz2,ntimes);
        nCCP = nan(nz2,ntimes);
        nDN = nan(nz2,ntimes);
        nTH = nan(nz2,ntimes);

        rams_hdf5c({'fthrd','vapliqt','nuccldrt','rcp','rrp','dn0',...
            'theta','wp','ccp'},times,folders{i,ff});

        for t=1:ntimes
            %find BL top
            th1h=squeeze(mean(mean(runs.theta(:,:,:,t))));
            th1h=gradient(th1h);
            [~,indz]=max(th1h);

            bottom = CTind - indz + 1;
            top = bottom - 1 + nz;
            nFTH(bottom:top,t) = squeeze(mean(mean(runs.fthrd(:,:,:,t))));
            nVAP(bottom:top,t) = squeeze(mean(mean(runs.vapliqt(:,:,:,t)+...
                runs.nuccldrt(:,:,:,t))))/(15*60);%convert to kg/kg/s
            nRCP(bottom:top,t) = squeeze(mean(mean(runs.rcp(:,:,:,t))));
            nCCP(bottom:top,t) = squeeze(mean(mean(runs.ccp(:,:,:,t))));
            nRTP(bottom:top,t) = squeeze(mean(mean(runs.rcp(:,:,:,t)+...
                runs.rrp(:,:,:,t))));
            nDN(bottom:top,t) = squeeze(mean(mean(runs.dn0(:,:,:,t))));
            nTH(bottom:top,t) = squeeze(mean(mean(runs.theta(:,:,:,t))));
        end

        mnDN=mean(nDN,2);
        mnFTH=mean(nFTH.*nDN,2);
        radcool(i,ff)=nansum(mnFTH(CTz<=0))*1004*dz;
        mnVAP=mean(nVAP.*nDN,2);

        for ic=1:length(ctd)
            evaptop(i,ff,ic)=sum(mnVAP(CTz>=ctd(ic)&CTz<=10))*2.5e6*dz;
            radtop(i,ff,ic)=sum(mnFTH(CTz>=ctd(ic)&CTz<=10))*1004*dz;
        end

        mLWP(i,ff)=nansum(mean(nRTP.*nDN,2))*dz;
        mFTH(:,i,ff)=mean(nFTH,2)*3600; %K/hr
        mRCP(:,i,ff)=mean(nRCP,2);
        mCCP(:,i,ff)=mean(nCCP,2);
        mRTP(:,i,ff)=mean(nRTP,2);
        mVAP(:,i,ff)=mean(nVAP,2)*2.5e6/1004*3600; %K/hr
        mTH(:,i,ff)=mean(nTH,2);
    end
end

save('isdacNew_dycomsNewSubs_hr_5.mat','we','we2','evaptop','radtop',...
    'radcool','mLWP','mFTH','mRCP','mVAP','mCCP','mRTP','mTH','entheat')

%% FIGURE 1
figure('position',[1350 300 1050 700])
t=tiledlayout(length(folds),2,'TileSpacing','tight');

load('isdacNew_dycomsNewSubs5.mat'); %low res sims
for cd=1:2
    for ff=1:length(folds)
        nexttile
        hold on
        for i=1:8
            for j=1:2
                h=plot((we(i,j,ff,cd))*100,mLWP(i,j,ff,cd)*1000,...
                    symsinv(i),'LineWidth',lwinv(i),...
                    'MarkerSize',10,'Color',cols((2-mod(ff,2))+2*(j-1),:),...
                    'MarkerFaceColor',cols((2-mod(ff,2))+2*(j-1),:))  ;
                if j==2 ; legendlabel(h);end
            end
        end
        grid on
        if ff==1
            set(gca,'xlim',[0.38 0.44],'ylim',[26 43])
            if cd==1
                title(['DYCOMS',10,'\rm\fontsize{15}                    ',...
                    '                                   ',...
                    'Tests Use "Clean" Diameter'])
                lg=legend(lgtextinv);
                lg.Layout.Tile='east';
                title(lg,[lgtitle,'Processes Modified:'])
            else
                title(['\rm\fontsize{15}                                ',...
                    '                         ',...
                    'Tests Use "Polluted" Diameter'])
            end

        elseif ff==2
            set(gca,'xlim',[0.17 0.3],'ylim',[25 35])
            if cd==1
                title(['ISDAC',10,''])
            end
        end
        plotlabel(['(',alphabet(ff+(cd-1)*2),')'],'upperright')
    end
end
xlabel(t,'Entrainment Rate (cm s^{-1})','FontSize',18)
ylabel(t,'Liquid Water Path (g m^{-2})','FontSize',18)

my_savefig('cloud_top_cooling_final/Figure1',gcf)
%% FIGURE 2
load isdacNew_dycomsNewSubs_hr_5.mat

CTind=450;
dz=3;
CTz=([1:900]-CTind)*dz;

ylm=[-100 15];
figure('position',[1000 275 1200 800])
tp=tiledlayout(length(folds),3,'TileSpacing','tight');

for ff=1:length(folds)
    nexttile;hold on
    for i=[1:6]
        if i==6
            plot(mVAP(:,i,ff),CTz,':','color',[.3 .4 .1])
        else
            plot(mVAP(:,i,ff),CTz,'-')
        end

        set(gca,'ylim',ylm)
        if ff==1
            set(gca,'xlim',[-7 1.5])
        else
            set(gca,'xlim',[-2.5 0.75])
        end
    end
    ylabel('Dist. from BL Top (m)')
    grid on
    if ff==2;xlabel('Lat. Heat. Rate (K hr^{-1})');end
    if ff==1
        h=legend({'None (DC\_M000)','Evap. Only (DC\_M100P)','Sed. Only (DC\_M010P)',...
            'Rad. Only (DC\_M001P)','All (DC\_M111P)','None (DP\_M000)'},...
            'Location','southwest','FontSize',12);
        title(h,'Processes Modified:')
    else
        h=legend({'None (IC\_M000)','Evap. Only (IC\_M100P)','Sed. Only (IC\_M010P)',...
            'Rad. Only (IC\_M001P)','All (IC\_M111P)','None (IP\_M000)'},...
            'Location','southwest','FontSize',12);
        title(h,'Processes Modified:')
    end
    plotlabel(['(',alphabet((ff-1)*3+1),')'])

    nexttile;hold on
    for i=[1:6]
        if i==6
            plot(mFTH(:,i,ff),CTz,':','color',[.3 .4 .1])
        else
            plot(mFTH(:,i,ff),CTz,'-')
        end

        set(gca,'ylim',ylm)
        if ff==1
            set(gca,'xlim',[-7 1.5])
        else
            set(gca,'xlim',[-2.5 0.75])
        end
    end
    set(gca,'yticklabel','')
    grid on
    if ff==2;xlabel('Rad. Heat. Rate (K hr^{-1})');end
    plotlabel(['(',alphabet((ff-1)*3+2),')'])

    nexttile;hold on
    for i=[1:6]
        if i==6
            plot(mRTP(:,i,ff)*1000,CTz,':','color',[.3 .4 .1])
        else
            plot(mRTP(:,i,ff)*1000,CTz,'-')
        end

        set(gca,'ylim',ylm)
        if ff==1
            set(gca,'xlim',[0 0.4])
        else
            set(gca,'xlim',[0 0.2])
        end
    end
    set(gca,'yticklabel','')
    grid on
    if ff==2;xlabel('Cloud Water (g kg^{-1})');end
    plotlabel(['(',alphabet((ff-1)*3+3),')'])
end
my_savefig('cloud_top_cooling_final/Figure2',gcf)
%% FIGURE 3&5

for itype=1:2
    load('isdacNew_dycomsNewSubs5.mat'); %low res sims
    figure('position',[1350 750 925 475])
    t=tiledlayout(length(folds),2,'TileSpacing','compact');
    jj=0;
    polors=[1 1 1;.75 0 0; .75 .75 0; 0 .75 0; 0 .75 .75; 0 0 .75; .75 0 .75;.75 0 0;0 0 0];
    for cd=1:2
        for ff=1:length(folds)
            jj=jj+1;
            axs(jj)=nexttile;
            hold on
            for i=1:8

                j=3-cd;
                pm=-1*(j*2-3);

                iord=[1,3,5,4,6,2,7,3,8];
                for ip=2:9
                    if itype==1
                        b1=bar(ip,pm*100*(we(iord(ip),j,ff,cd)-we(1,j,ff,cd)));
                        if any([5:7]==iord(ip))
                            h=plot(ip,pm*100*(we(iord(ip-1),j,ff,cd)-we(1,j,ff,cd)+...
                                we(iord(ip+1),j,ff,cd)-we(1,j,ff,cd)),...
                                'k*','MarkerSize',10,'LineWidth',1);
                            legendlabel(h)
                        end
                    else
                        b1=bar(ip,pm*1000*(mLWP(iord(ip),j,ff,cd)-mLWP(1,j,ff,cd)));
                        if any([5:7]==iord(ip))
                            h=plot(ip,pm*1000*(mLWP(iord(ip-1),j,ff,cd)-mLWP(1,j,ff,cd)+...
                                mLWP(iord(ip+1),j,ff,cd)-mLWP(1,j,ff,cd)),...
                                'k*','MarkerSize',10,'LineWidth',1);
                            legendlabel(h)
                        end
                    end
                    hold on
                    if ip>1
                        b1.FaceColor='flat';
                        b1.CData=polors(ip,:);
                    end
                end
            end
            grid on
            set(gca,'xtick',[])
            if itype==1
                plotlabel(['(',alphabet(jj),')'])
            else
                plotlabel(['(',alphabet(jj),')'],'lowerleft')
            end
            if ff==1
                if cd==1
                    title(['DYCOMS',10,'\rm\fontsize{15}                ',...
                        '                                   ',...
                        'Tests Use "Clean" Diameter'])
                end
                if cd==2
                    title(['\rm\fontsize{15}                           ',...
                        '                         ',...
                        'Tests Use "Polluted" Diameter'])
                end
            elseif ff==2
                if cd==1; title(['ISDAC',10]);end
            end

            if itype==1
                ylabel(t,'\Deltaw_e (cm s^{-1})','FontSize',18)
                if ff==1
                    set(gca,'ylim',[-0.005 0.05])
                else
                    set(gca,'ylim',[-0.02 0.12])
                end
            else
                ylabel(t,'\DeltaLWP (g m^{-2})','FontSize',18)
                if ff==1
                    set(gca,'ylim',[-20 0],'ytick',-20:5:0)
                else
                    set(gca,'ylim',[-12 0],'ytick',[-12:3:0])
                end
            end
        end
    end

    load('isdacNew_dycomsNewSubs_hr_5.mat'); %high res sims
    j=0;
    polors=[1 1 1;1 0 0;  0 1 0;  0 0 1; 1 0 0; .25 .25 .25];
    for cd=1:2
        for ff=1:length(folds)
            j=j+1;
            axes(axs(j))
            for i=1:5


                pm=cd*2-3;
                iord=(3-cd)*5-[5,3,2,4,3,1]+1;
                for ip=2:6
                    ipx=(ip-1)*2;
                    if ip==6; ipx=9;end
                    if itype==1
                        b1=bar(ipx,pm*100*(we(iord(ip),ff)-we(iord(1),ff)),'BarWidth',0.4);
                    else
                        b1=bar(ipx,pm*1000*(mLWP(iord(ip),ff)-mLWP(iord(1),ff)),'BarWidth',0.4);
                    end
                    hold on
                    if ip>1
                        b1.FaceColor='flat';
                        b1.CData=polors(ip,:);
                    end
                end
            end
        end
    end
    lg=legend(lgtextpolarinv{2:end});
    lg.Layout.Tile='east';
    title(lg,['Processes',10,'Changed'])

    if itype==1
        my_savefig('cloud_top_cooling_final/Figure3',gcf)
    else
        my_savefig('cloud_top_cooling_final/Figure5',gcf)
    end
end


%% FIGURE 4

figure('position',[1000 700 1150 900])

icl=[3,6];
ich=[3,6];

for ff=1:length(folds)

    axes('position',apos(ff,:))
    hold on
    xlim padded
    ylim padded
    load('isdacNew_dycomsNewSubs5.mat'); %low res sims
    ic=icl(ff);
    for j=1:2
        for i=[2:4,8]
            pm=(j*2-3);

            h=plot(pm*(radtop(i,j,ff,3-j,ic)-radtop(1,j,ff,3-j,ic))+...
                pm*(evaptop(i,j,ff,3-j,ic)-evaptop(1,j,ff,3-j,ic)),...
                -1*pm*100*(we(i,j,ff,3-j)-we(1,j,ff,3-j)),symsinv(i),'LineWidth',lwinv(i),...
                'MarkerSize',10,'Color',cols((2-mod(ff,2))+2*(j-1),:),...
                'MarkerFaceColor',cols((2-mod(ff,2))+2*(j-1),:));
        end
    end
    if ff==1
        title('DYCOMS')
        lg=legend(lgtextdiff);
        lg.Position=[.7733 .2 .2174 .55];
        title(lg,[lgtitle2,'Processes Modified:'])
        plotlabel('(a)')
        xticks(-2:2:10)
    end

    load('isdacNew_dycomsNewSubs_hr_5.mat'); %low res sims
    ic=ich(ff);
    for i=[2:5,7:10]
        isty=i-floor((i-1)/5)*5;
        pm=(floor((i-1)/5)*2-1);
        ib=floor((i-1)/5)*5+1;
        h=plot(pm*(radtop(i,ff,ic)-radtop(ib,ff,ic))+pm*(evaptop(i,ff,ic)-evaptop(ib,ff,ic)),...
            -1*pm*100*(we(i,ff)-we(ib,ff)),symsinvhr(isty),'LineWidth',lwinvhr(isty)^2/1.5,...
            'MarkerSize',14,'Color','k','MarkerFaceColor','k');
        legendlabel(h)
        h=plot(pm*(radtop(i,ff,ic)-radtop(ib,ff,ic))+pm*(evaptop(i,ff,ic)-evaptop(ib,ff,ic)),...
            -1*pm*100*(we(i,ff)-we(ib,ff)),symsinvhr(isty),'LineWidth',lwinvhr(isty),...
            'MarkerSize',8,'Color',cols((2-mod(ff,2))+floor((i-1)/5)*2,:),...
            'MarkerFaceColor',cols((2-mod(ff,2))+floor((i-1)/5)*2,:));
        legendlabel(h)
    end

    grid on
    h=plot(get(gca,'xlim'),[0 0],'--','LineWidth',0.5,'Color',[.5 .5 .5]);
    legendlabel(h)
    h=plot([0 0],get(gca,'ylim'),'--','LineWidth',0.5,'Color',[.5 .5 .5]);
    legendlabel(h)

    xlabel('\Delta CT Diab. Cooling (W m^{-2})')
    if ff==1;ylabel('\Delta w_e (cm s^{-1})');end
    if ff==2;title('ISDAC');plotlabel('(b)');xticks(-1:10);end
end
for ff=1:length(folds)

    axes('position',apos(ff+2,:))
    hold on
    xlim padded
    ylim padded
    load('isdacNew_dycomsNewSubs5.mat'); %low res sims
    ic=icl(ff);
    for i=[2:4,8]
        for j=1:2
            pm=(j*2-3);
            h=plot(pm*(evaptop(i,j,ff,3-j,ic)-evaptop(1,j,ff,3-j,ic)),...
                pm*(radtop(i,j,ff,3-j,ic)-radtop(1,j,ff,3-j,ic)),...
                symsinv(i),'LineWidth',lwinv(i),...
                'MarkerSize',10,'Color',cols((2-mod(ff,2))+2*(j-1),:),...
                'MarkerFaceColor',cols((2-mod(ff,2))+2*(j-1),:));

            if j==2 ; legendlabel(h);end
        end
    end
    if ff==1
        plotlabel('(c)')
    end
    load('isdacNew_dycomsNewSubs_hr_5.mat'); %low res sims
    ic=ich(ff);
    for i=[2:5,7:10]
        isty=i-floor((i-1)/5)*5;
        pm=(floor((i-1)/5)*2-1);
        ib=floor((i-1)/5)*5+1;
        h=plot(pm*(evaptop(i,ff,ic)-evaptop(ib,ff,ic)),...
            pm*(radtop(i,ff,ic)-radtop(ib,ff,ic)),symsinvhr(isty),...
            'LineWidth',lwinvhr(isty)^2/1.5,...
            'MarkerSize',14,'Color','k','MarkerFaceColor','k');
        legendlabel(h)
        h=plot(pm*(evaptop(i,ff,ic)-evaptop(ib,ff,ic)),...
            pm*(radtop(i,ff,ic)-radtop(ib,ff,ic)),symsinvhr(isty),...
            'LineWidth',lwinvhr(isty),...
            'MarkerSize',8,'Color',cols((2-mod(ff,2))+floor((i-1)/5)*2,:),...
            'MarkerFaceColor',cols((2-mod(ff,2))+floor((i-1)/5)*2,:));
        legendlabel(h)
    end

    grid on
    h=plot(get(gca,'xlim'),[0 0],'--','LineWidth',0.5,'Color',[.5 .5 .5]);
    legendlabel(h)
    h=plot([0 0],get(gca,'ylim'),'--','LineWidth',0.5,'Color',[.5 .5 .5]);
    legendlabel(h)

    xlabel('\Delta CT Evap. Cooling (W m^{-2})')
    if ff==1;ylabel('\Delta CT Rad. Cooling (W m^{-2})');end
    if ff==2;plotlabel('(d)');yticks(-1:4);end
end
my_savefig('cloud_top_cooling_final/Figure4',gcf)


%% FIGURE 6

figure('position',[1000 700 1150 900])

for ff=1:length(folds)

    axes('position',apos(ff,:))
    hold on
    xlim padded
    ylim padded
    load('isdacNew_dycomsNewSubs5.mat'); %low res sims

    for j=1:2
        for i=[2:4,8]
            pm=-1*(j*2-3);
            h=plot(pm*(entheat(i,j,ff,3-j)-entheat(1,j,ff,3-j))+pm*(radcool(i,j,ff,3-j)-radcool(1,j,ff,3-j)),...
                pm*1000*(mLWP(i,j,ff,3-j)-mLWP(1,j,ff,3-j)),symsinv(i),'LineWidth',lwinv(i),...
                'MarkerSize',10,'Color',cols((2-mod(ff,2))+2*(j-1),:),'MarkerFaceColor',cols((2-mod(ff,2))+2*(j-1),:));
        end
    end

    if ff==1
        title('DYCOMS')
        lg=legend(lgtextdiff);
        lg.Position=[.7733 .2 .2174 .55];
        title(lg,[lgtitle2,'Processes Modified:'])
        plotlabel('(a)','upperright')
    end
    load('isdacNew_dycomsNewSubs_hr_5.mat'); %low res sims

    for i=[2:5,7:10]
        isty=i-floor((i-1)/5)*5;
        pm=-1*(floor((i-1)/5)*2-1);
        ib=floor((i-1)/5)*5+1;
        h=plot(pm*(entheat(i,ff)-entheat(ib,ff))+pm*(radcool(i,ff)-radcool(ib,ff)),...
            pm*1000*(mLWP(i,ff)-mLWP(ib,ff)),symsinvhr(isty),'LineWidth',lwinvhr(isty)^2/1.5,...
            'MarkerSize',14,'Color','k','MarkerFaceColor','k');
        legendlabel(h)
        h=plot(pm*(entheat(i,ff)-entheat(ib,ff))+pm*(radcool(i,ff)-radcool(ib,ff)),...
            pm*1000*(mLWP(i,ff)-mLWP(ib,ff)),symsinvhr(isty),'LineWidth',lwinvhr(isty),...
            'MarkerSize',8,'Color',cols((2-mod(ff,2))+floor((i-1)/5)*2,:),...
            'MarkerFaceColor',cols((2-mod(ff,2))+floor((i-1)/5)*2,:));
        legendlabel(h)
    end

    grid on
    h=plot(get(gca,'xlim'),[0 0],'--','LineWidth',0.5,'Color',[.5 .5 .5]);
    legendlabel(h)
    h=plot([0 0],get(gca,'ylim'),'--','LineWidth',0.5,'Color',[.5 .5 .5]);
    legendlabel(h)

    xlabel('\Delta(H_{ent}+H_{rad}) (W m^{-2})')
    if ff==1;ylabel('\Delta LWP (g kg^{-1})');end
    if ff==2;title('ISDAC');plotlabel('(b)','upperright');end
end
for ff=1:length(folds)

    axes('position',apos(ff+2,:))
    hold on
    xlim padded
    ylim padded
    load('isdacNew_dycomsNewSubs5.mat'); %low res sims

    for i=[2:4,8]
        for j=1:2
            pm=-1*(j*2-3);

            h=plot(pm*(entheat(i,j,ff,3-j)-entheat(1,j,ff,3-j)),...
                -1*pm*(radcool(i,j,ff,3-j)-radcool(1,j,ff,3-j)),...
                symsinv(i),'LineWidth',lwinv(i),...
                'MarkerSize',10,'Color',cols((2-mod(ff,2))+2*(j-1),:),...
                'MarkerFaceColor',cols((2-mod(ff,2))+2*(j-1),:));
            if j==2 ; legendlabel(h);end
        end
    end

    if ff==1
        plotlabel('(c)')
    end
    load('isdacNew_dycomsNewSubs_hr_5.mat'); %low res sims

    for i=[2:5,7:10]
        isty=i-floor((i-1)/5)*5;
        pm=-1*(floor((i-1)/5)*2-1);
        ib=floor((i-1)/5)*5+1;
        h=plot(pm*(entheat(i,ff)-entheat(ib,ff)),-1*pm*(radcool(i,ff)-...
            radcool(ib,ff)),symsinvhr(isty),'LineWidth',lwinvhr(isty)^2/1.5,...
            'MarkerSize',14,'Color','k','MarkerFaceColor','k');
        legendlabel(h)
        h=plot(pm*(entheat(i,ff)-entheat(ib,ff)),-1*pm*(radcool(i,ff)-...
            radcool(ib,ff)),symsinvhr(isty),'LineWidth',lwinvhr(isty),...
            'MarkerSize',8,'Color',cols((2-mod(ff,2))+floor((i-1)/5)*2,:),...
            'MarkerFaceColor',cols((2-mod(ff,2))+floor((i-1)/5)*2,:));
        legendlabel(h)
    end

    grid on
    h=plot(get(gca,'xlim'),[0 0],'--','LineWidth',0.5,'Color',[.5 .5 .5]);
    legendlabel(h)
    h=plot([0 0],get(gca,'ylim'),'--','LineWidth',0.5,'Color',[.5 .5 .5]);
    legendlabel(h)

    xlabel('\Delta Entrainment Heating (W m^{-2})')
    if ff==1;ylabel('\Delta Radiative Cooling (W m^{-2})');end
    if ff==2;plotlabel('(d)');end
end
my_savefig('cloud_top_cooling_final/Figure6',gcf)

%% FIGURE 7

figure('position',[1000 700 1250 500])
t=tiledlayout(1,2,'TileSpacing','compact');
for ff=1:length(folds)

    nexttile
    hold on
    xlim padded
    ylim padded
    load('isdacNew_dycomsNewSubs5.mat'); %low res sims

    for j=1:2
        for i=[2:4,8]
            h=plot(1000*mLWP(i,j,ff,3-j),-1*radcool(i,j,ff,3-j),symsinv(i),'LineWidth',lwinv(i),...
                'MarkerSize',10,'Color',cols((2-mod(ff,2))+2*(j-1),:),...
                'MarkerFaceColor',cols((2-mod(ff,2))+2*(j-1),:));
        end
    end

    if ff==1
        title('DYCOMS')
        lg=legend(lgtextdiff);
        lg.Layout.Tile='east';
        title(lg,[lgtitle2,'Processes Modified:'])
        plotlabel('(a)')
    end
    load('isdacNew_dycomsNewSubs_hr_5.mat'); %low res sims

    for i=[2:5,7:10]
        isty=i-floor((i-1)/5)*5;
        pm=-1*(floor((i-1)/5)*2-1);
        ib=floor((i-1)/5)*5+1;
        h=plot(1000*mLWP(i,ff),-1*radcool(i,ff),symsinvhr(isty),'LineWidth',lwinvhr(isty)^2/1.5,...
            'MarkerSize',14,'Color','k','MarkerFaceColor','k');
        legendlabel(h)
        h=plot(1000*mLWP(i,ff),-1*radcool(i,ff),symsinvhr(isty),'LineWidth',lwinvhr(isty),...
            'MarkerSize',8,'Color',cols((2-mod(ff,2))+floor((i-1)/5)*2,:),...
            'MarkerFaceColor',cols((2-mod(ff,2))+floor((i-1)/5)*2,:));
        legendlabel(h)
    end

    grid on

    xlabel('LWP (g kg^{-1})')
    if ff==1;ylabel('Radiative Cooling (W m^{-2})');end
    if ff==2;title('ISDAC');plotlabel('(b)');end
end

my_savefig('cloud_top_cooling_final/Figure7',gcf)
