folds={'dycoms','isdac'};
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)']};
lgtextdiffsimple={['Evap. Only - C\_M000'],['Sed. Only - C\_M000'],...
    ['Rad. Only - C\_M000'],['All - C\_M000'],...
    ['P\_M000 - Evap. Only'],['P\_M000 - Sed. Only'],...
    ['P\_M000 - Rad. Only'],['P\_M000 - All']};
lgtextdiffsimple2={['Evap. Only/C\_M000'],['Sed. Only/C\_M000'],...
    ['Rad. Only/C\_M000'],['All/C\_M000'],...
    ['P\_M000/Evap. Only'],['P\_M000/Sed. Only'],...
    ['P\_M000/Rad. Only'],['P\_M000/All']};
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,diabheat,diabheat2,mass,mLWP,entheat,entwater,effic,effic2,sflx,lflx,wstar,ziavg]=deal(nan(8,2,length(folds),2));
[LWPts,zits]=deal(nan(73,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));
load('isdac_dycoms.mat')
for cd=1:2 %fixed diameter
    for ff=1:length(folds) %dycoms, isdac
        if ff==1
            times=12:48;
        else
            times=8:72;
        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},'_000_0']; %full physics
        folders{1,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_000_0'];
        folders{2,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_020_',num2str(cdiam)]; %noevap
        folders{2,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_020_',num2str(cdiam)];
        folders{3,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_100_',num2str(cdiam)]; %nosed
        folders{3,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_100_',num2str(cdiam)];
        folders{4,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_001_',num2str(cdiam)]; %norad
        folders{4,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_001_',num2str(cdiam)];
        folders{5,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_101_',num2str(cdiam)]; %evap only
        folders{5,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_101_',num2str(cdiam)];
        folders{6,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_021_',num2str(cdiam)]; %sed only
        folders{6,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_021_',num2str(cdiam)];
        folders{7,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_120_',num2str(cdiam)]; %rad only
        folders{7,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_120_',num2str(cdiam)];
        folders{8,1,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_121_',num2str(cdiam)]; %deny all
        folders{8,2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_121_',num2str(cdiam)];

        for i=1:8 %physics tests
            for j=1:2 %aerosol conc
                display([cd,ff,i,j])

                rams_hdf5c({'theta','rv','rcp','pi','dn0','rtp'},[0:times(end)],folders{i,j,ff});

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

                for t=1:size(runs.tempk,4)
                    th1=squeeze(mean(mean(runs.theta(:,:,:,t))));

                    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

                    zits(t,i,j,ff,cd)=z1i; %BL height change - average subsidence rate

                end
                entrainment rate calculations - doing two ways
                clear global

                rams_hdf5c({'theta','rv','rcp','pi','dn0','rtp'},[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
                indzkeep=indz;
                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; %seconds between start and end
                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 heat and water fluxes
                bltind=round(indzkeep/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;

                meanQ=squeeze(mean(mean(runs.rtp(:,:,:,1))));
                blQ=trapz(meanQ(1:bltind).*squeeze(runs.dn0(1,1,1:bltind,1)))/trapz(runs.dn0(1,1,1:bltind,1));
                dQ=min(meanQ)-blQ;

                dntop=runs.dn0(1,1,bltind);

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

                ziavg(i,j,ff,cd)=(z2i+z1i)/2;
                rams_hdf5c({'wp_buoy_theta','wp_buoy_cond','wp','sflux_t','sflux_r'},times,folders{i,j,ff});
                bp=runs.wp_buoy_theta+runs.wp_buoy_cond;
                zim=(z2i+z1i)/2;
                izim=round((indz+indzkeep)/200);
                dbp = mean(min(min(min(bp))) - mean(mean(mean(bp(:,:,1:bltind,:)))));
                bp=bp-mean(mean(bp));

                for t=1:size(runs.wp,4)
                    wpbp=squeeze(mean(mean(bp(:,:,:,t).*runs.wp(:,:,:,t))));
                    wst(t)=mean(2.5*(trapz(wpbp)*dz)).^(1/3);
                end
                effic(i,j,ff,cd)=-1*we(i,j,ff,cd)*zim/mean(wst.^3)*dbp;
                effic2(i,j,ff,cd)=-1*we(i,j,ff,cd)*zim/mean(wst).^3*dbp;
                wstar(i,j,ff,cd)=mean(wst);
                sflx(i,j,ff,cd)=mean(mean(mean(runs.sflux_t)))*1004;
                lflx(i,j,ff,cd)=mean(mean(mean(runs.sflux_r)))*2.5e6;

                %BL-top normalized profiles and integrals
                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);

                clear global
                rams_hdf5c({'fthrd','vapliqt','nuccldrt','rcp','rrp','dn0',...
                    'theta','wp','ccp'},times,folders{i,j,ff});
 
                nz=size(runs.theta,3);
                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,j,ff,cd)=nansum(mnFTH(CTz<=0))*1004*dz;

                mnVAP=mean(nVAP.*nDN,2);
                diabheat(i,j,ff,cd)=nansum(mnFTH(CTz<=0)*1004*dz+mnVAP(CTz<=0)*2.5e6*dz);
                diabheat2(i,j,ff,cd)=nansum(mnFTH(CTz<=10)*1004*dz+mnVAP(CTz<=10)*2.5e6*dz);

                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);

                clear global
                rams_hdf5c({'rcp','rrp','dn0'},0:times(end),folders{i,j,ff});
                ntimes2=size(runs.rcp,4);
                LWPts(1:ntimes2,i,j,ff,cd)=squeeze(sum(mean(mean((runs.rcp+runs.rrp).*runs.dn0)),3)*dz);

            end
        end
    end
end
save('isdac_dycoms.mat','we','we2','evaptop','radtop','radcool','LWPts','diabheat','diabheat2',...
    'mLWP','mFTH','mRCP','mVAP','mCCP','mRTP','mTH','entheat','entwater','effic','zits',...
    'effic2','sflx','lflx','wstar','ziavg')

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

%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,entwater,effic,effic2,sflx,lflx,wstar,ziavg]=deal(nan(10,length(folds)));
[LWPts,zits]=deal(nan(73,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));
load isdac_dycoms_highres
for ff=2%1:length(folds)
    if ff==1
        times=12:48;
    else
        times=8: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},'_000_0']; %full physics
    folders{2,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_020_',num2str(cdiams(ff,2))]; %noevap
    folders{3,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_100_',num2str(cdiams(ff,2))]; %nosed
    folders{4,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_001_',num2str(cdiams(ff,2))]; %norad
    folders{5,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_',ccn_clean{ff},'_121_',num2str(cdiams(ff,2))]; %deny all
    folders{6,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_000_0'];
    folders{7,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_020_',num2str(cdiams(ff,1))]; %noevap
    folders{8,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_100_',num2str(cdiams(ff,1))]; %nosed
    folders{9,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_001_',num2str(cdiams(ff,1))]; %norad
    folders{10,ff}=['~/rams-Arctic-Abbey/',folds{ff},'/output_1600_121_',num2str(cdiams(ff,1))]; %deny all

    for i=1:10
        toc
        tic
        [ff,i]
        rams_hdf5c({'theta','rv','rcp','pi','dn0','rtp'},[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
        indzkeep=indz;
        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);

        ziavg(i,ff)=(z2i+z1i)/2;

        bltind=round(indzkeep/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;

        meanQ=squeeze(mean(mean(runs.rtp(:,:,:,1))));
        blQ=trapz(meanQ(1:bltind).*squeeze(runs.dn0(1,1,1:bltind,1)))/trapz(runs.dn0(1,1,1:bltind,1));
        dQ=min(meanQ)-blQ;

        dntop=runs.dn0(1,1,bltind);

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

        rams_hdf5c({'wp_buoy_theta','wp_buoy_cond','wp','sflux_t','sflux_r'},times,folders{i,ff});
        bp=runs.wp_buoy_theta+runs.wp_buoy_cond;
        zim=(z2i+z1i)/2;
        izim=round((indz+indzkeep)/200);
        dbp = mean(mean(mean(mean(bp(:,:,1:izim,:)))))-mean(mean(mean(bp(:,:,izim+1,:))));
        bp=bp-mean(mean(bp));

        for t=1:size(runs.wp,4)
            wpbp=squeeze(mean(mean(bp(:,:,:,t).*runs.wp(:,:,:,t))));
            wst(t)=mean(2.5*(trapz(wpbp)*dz)).^(1/3);
        end

        effic(i,ff)=we(i,ff)*zim/mean(wst.^3)*dbp;
        effic2(i,ff)=we(i,ff)*zim/mean(wst).^3*dbp;
        wstar(i,ff)=mean(wst);
        sflx(i,ff)=mean(mean(mean(runs.sflux_t)))*1004;
        lflx(i,ff)=mean(mean(mean(runs.sflux_r)))*2.5e6;

        %BL-top normalized profiles and integrals


        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);

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

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

            bottom = CTind - indz + 1;
            top = bottom - 1 + nz;

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

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

        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);
        clear global
        rams_hdf5c({'rcp','rrp','dn0'},0:times(end),folders{i,ff});
        ntimes2=size(runs.rcp,4);
        LWPts(1:ntimes2,i,ff)=squeeze(sum(mean(mean((runs.rcp+runs.rrp).*runs.dn0)),3)*dz);
    end
end

save('isdac_dycoms_highres.mat','we','we2','evaptop','radtop','effic','effic2','wstar','sflx','lflx','diabheat','diabheat2',...
    'radcool','mLWP','mFTH','mRCP','mVAP','mCCP','mRTP','mTH','entheat','entwater','ziavg','LWPts')

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

load('isdac_dycoms.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
            if cd==1
                title(['DYCOMS',10,'\rm\fontsize{15}                    ',...
                    '                             ',...
                    'Tests Use "Clean" Diameter (MxxxC Tests)'])
                lg=legend(lgtextinv);
                lg.Layout.Tile='east';
                title(lg,[lgtitle,'Processes Modified:'])
            else
                title(['\rm\fontsize{15}                                ',...
                    '                   ',...
                    'Tests Use "Polluted" Diameter (MxxxP Tests)'])
            end

        elseif ff==2
            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 isdac_dycoms_highres.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&6

for itype=1:2
    load('isdac_dycoms.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;
            if ff==1;iend=48;else;iend=60;end
            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)),...
                                'o','MarkerSize',10,'LineWidth',1,'MarkerFaceColor',polors(ip,:),'MarkerEdgeColor','w');
                            legendlabel(h)
                        elseif iord(ip)==8
                            wesave(cd,ff)=pm*100*sum(we(2:4,j,ff,cd)-we(1,j,ff,cd));
                        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)),...
                                'o','MarkerSize',10,'LineWidth',1,'MarkerFaceColor',polors(ip,:),'MarkerEdgeColor','w');
                            legendlabel(h)
                        elseif iord(ip)==8
                            %don't plot this now, otherwise it'll be
                            %covered by high-res bars
                            lwpsave(cd,ff)=pm*1000*sum(mLWP(2:4,j,ff,cd)-mLWP(1,j,ff,cd));
                        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 .05])
                else
                    set(gca,'ylim',[-0.01 0.12])
                end
            else
                ylabel(t,'\DeltaLWP (g m^{-2})','FontSize',18)
                if ff==1
                    set(gca,'ylim',[-25 0],'ytick',-25:5:0)
                else
                    set(gca,'ylim',[-12 0],'ytick',[-12:3:0])
                end
            end
        end
    end

    load('isdac_dycoms_highres.mat'); %high res sims

    j=0;
    polors=[1 1 1;1 0 0;  0 1 0;  0 0 1; 1 0 0; .5 .5 .5];
    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);
                        if ip==6
                            %plot sum of low-res simulations now
                            h=plot(ipx,wesave(cd,ff),...
                                'o','MarkerSize',10,'LineWidth',1,'MarkerEdgeColor','w','MarkerFaceColor','k');
                            legendlabel(h)

                            h=plot(ipx,pm*100*sum(we(iord(2:4),ff)-we(iord(1),ff)),...
                                'o','MarkerSize',7,'LineWidth',1,'MarkerFaceColor',polors(6,:),'MarkerEdgeColor','w');
                            legendlabel(h)
                        end
                    else
                        b1=bar(ipx,pm*1000*(mLWP(iord(ip),ff)-mLWP(iord(1),ff)),'BarWidth',0.4);

                        if ip==6
                            %plot sum of low-res simulations now
                            h=plot(ipx,lwpsave(cd,ff),...
                                'o','MarkerSize',10,'LineWidth',1,'MarkerFaceColor','k','MarkerEdgeColor','w');
                            legendlabel(h)

                            h=plot(ipx,pm*1000*sum(mLWP(iord(2:4),ff)-mLWP(iord(1),ff)),...
                                'o','MarkerSize',7,'LineWidth',1,'MarkerFaceColor',polors(6,:),'MarkerEdgeColor','w');
                            legendlabel(h)
                        end
                    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/Figure6',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('isdac_dycoms.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('isdac_dycoms_highres.mat'); %high 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('isdac_dycoms.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('isdac_dycoms_highres.mat'); %high 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 5

figure('position',[1350 750 1000 350])
t=tiledlayout(1,length(folds),'TileSpacing','compact');
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 ff=1:length(folds)
    axs(ff)=nexttile;
    hold on
    axis padded

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

        load('isdac_dycoms.mat')

        for i=[2:4,8]
            dA=(effic(i,j,ff,cd)/effic(1,j,ff,cd));
            dB=(wstar(i,j,ff,cd)/wstar(1,j,ff,cd))^3;

            b1=plot((dA)^pm,(dB)^pm,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

        load('isdac_dycoms_highres.mat')

        for i=[2:5,7:10]
            isty=i-floor((i-1)/5)*5;
            pm=(floor((i-1)/5)*2-1)*-1;
            ib=floor((i-1)/5)*5+1;

            dA=(effic(i,ff)/effic(ib,ff));
            dB=(wstar(i,ff)/wstar(ib,ff))^3;

            h=plot((dA)^pm,(dB)^pm,symsinvhr(isty),...
                'LineWidth',lwinvhr(isty)^2/1.5,...
                'MarkerSize',14,'Color','k','MarkerFaceColor','k');
            legendlabel(h)
            h=plot((dA)^pm,(dB)^pm,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
        set(gca,'xscale','log','yscale','log')
        xlabel('Efficiency Ratio')
        if ff==1;ylabel('w_*^3 Ratio');end
        set(gca,'xlim',[min([get(gca,'xlim'),1]) max([1,get(gca,'xlim')])],'ylim',[min([1,get(gca,'ylim')]) max([1,get(gca,'ylim')])])
        if ff==1
            [~,h]=contour(.8:.01:2.2,.5:.01:1.5,[[.8:.01:2.2]'*[.5:.01:1.5]]',.1:.05:2,'Color',[.8 .8 .8]);
            legendlabel(h)
            [~,h]=contour(.8:.01:2.2,.5:.01:1.5,[[.8:.01:2.2]'*[.5:.01:1.5]]',[1 1],'Color','k');
            legendlabel(h)

            lg=legend(lgtextdiffsimple2);
            lg.Layout.Tile='east';
            title('DYCOMS')
        else
            contour(.8:.01:2.5,.5:.01:1.5,[[.8:.01:2.5]'*[.5:.01:1.5]]',.1:.1:2,'Color',[.8 .8 .8])
            contour(.8:.01:2.5,.5:.01:1.5,[[.8:.01:2.5]'*[.5:.01:1.5]]',[1 1],'Color','k')
            title('ISDAC')
        end

    end
end
my_savefig('cloud_top_cooling_final/Figure5',gcf)

%% FIGURE 7

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

for ff=1:length(folds)

    axes('position',apos(ff,:))
    hold on
    xlim padded
    ylim padded
    load('isdac_dycoms.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('isdac_dycoms_highres.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 m^{-2})');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('isdac_dycoms.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('isdac_dycoms_highres.mat');  %high 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/Figure7',gcf)

%% FIGURE 8

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('isdac_dycoms.mat'); %low res sims

    for j=1:2
        for i=[1: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(lgtextinv{[1:4,8]});
        lg.Layout.Tile='east';
        title(lg,[lgtitle2,'Processes Modified:'])
        plotlabel('(a)')
    end
    load('isdac_dycoms_highres.mat'); %low res sims

    for i=[1:5,6: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/Figure8',gcf)
%% FIGURE 9

figure('position',[1000 700 700 400])

ff=1;

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

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

        h=plot(pm*(sflx(i,j,ff,3-j)-sflx(1,j,ff,3-j)),...
            pm*(lflx(i,j,ff,3-j)-lflx(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

lg=legend(lgtextdiffsimple,'Location','eastoutside' );

load('isdac_dycoms_highres.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*(sflx(i,ff)-sflx(ib,ff)),...
        pm*(lflx(i,ff)-lflx(ib,ff)),symsinvhr(isty),'LineWidth',lwinvhr(isty)^2/1.5,...
        'MarkerSize',14,'Color','k','MarkerFaceColor','k');
    legendlabel(h)
    h=plot(pm*(sflx(i,ff)-sflx(ib,ff)),...
        pm*(lflx(i,ff)-lflx(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

xlim([-2 ,.25])
ylim([-.25,3.5])
set(gca,'ytick',0:3,'xtick',-2:0)
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 Sensible Heat Flux (W m^{-2})')
ylabel('\Delta Latent Heat Flux (W m^{-2})')

my_savefig('cloud_top_cooling_final/Figure9',gcf)
%%}