% % clear all; close all;
% % 
% % fil1 = 'C:\Users\Aaron Donohoe\Desktop\research\radiation\ceres_ebaf\4.1\CERES_EBAF_Ed4.1_Subset_200003-201608.nc';
% % fil2 = 'C:\Users\Aaron Donohoe\Desktop\research\radiation\ceres_ebaf\4.1\CERES_EBAF_Ed4.1_Subset_201609-202106.nc';
% % 
% % %ncdisp(fil1);
% % lon=double(ncread(fil1, 'lon',[1], [Inf])); N=length(lon);
% % lat=double(ncread(fil1, 'lat', [1], [Inf])); NN=length(lat);
% % t1=double(ncread(fil1, 'time', [1], [Inf])); NT1=length(t1); %days since 03 /01/2000
% % daysperyear = 365.2425;
% % 
% % t_off=monthday(3); %days to add to get in calandar days
% % y_off=2000;
% % year1=y_off+rounddown((t1+t_off)/daysperyear);
% % month1=roundup( ((t1+t_off) - daysperyear*(year1-y_off))/(daysperyear/12));
% % 
% % t2=double(ncread(fil2, 'time', [1], [Inf])); NT2=length(t2); %days since 03 /01/2000
% % 
% % year2=y_off+rounddown((t2+t_off)/daysperyear);
% % month2=roundup( ((t2+t_off) - daysperyear*(year2-y_off))/(daysperyear/12));
% % 
% % month=[month1; month2];
% % year=[year1; year2];
% % NT=NT1+NT2;
% % 
% % 
% % UT(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'toa_sw_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
% % UT((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'toa_sw_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
% %  
% % 
% % DT(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'solar_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
% % DT((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'solar_mon', [1 1 1], [Inf Inf Inf]), [3 2 1])); 
% % 
% % OLR(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'toa_lw_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
% % OLR((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'toa_lw_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1])); 
% % 
% % US(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'sfc_sw_up_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
% % US((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'sfc_sw_up_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1])); 
% % 
% % DS(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'sfc_sw_down_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
% % DS((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'sfc_sw_down_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1])); 
% % 
% % LDS(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'sfc_lw_down_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
% % LDS((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'sfc_lw_down_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1])); 
% % 
% % LUS(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'sfc_lw_up_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
% % LUS((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'sfc_lw_up_all_mon', [1 1 1], [Inf Inf Inf]), [3 2 1])); 
% % 
% % %make climatology and anomalies
% % for mn=1:12;
% %     do=find(month == mn);
% %     CDT(mn,1:NN,1:N) = squeeze(mean(DT(do,:,:)));
% %     CUT(mn,1:NN,1:N) = squeeze(mean(UT(do,:,:)));
% %     COLR(mn,1:NN,1:N) = squeeze(mean(OLR(do,:,:)));
% % 
% %     CDS(mn,1:NN,1:N) = squeeze(mean(DS(do,:,:)));
% %     CUS(mn,1:NN,1:N) = squeeze(mean(US(do,:,:)));
% %     CLDS(mn,1:NN,1:N) = squeeze(mean(LDS(do,:,:)));
% %     CLUS(mn,1:NN,1:N) = squeeze(mean(LUS(do,:,:)));
% % 
% %     for n=1:length(do);
% %      aDT(do(n), 1:NN,1:N) = DT(do(n),:,:)-CDT(mn,:,:);
% %      aUT(do(n), 1:NN,1:N) = UT(do(n),:,:)-CUT(mn,:,:);
% %      aOLR(do(n), 1:NN,1:N) = OLR(do(n),:,:)-COLR(mn,:,:);
% % 
% %      aDS(do(n), 1:NN,1:N) = DT(do(n),:,:)-CDT(mn,:,:);
% %      aUS(do(n), 1:NN,1:N) = US(do(n),:,:)-CUS(mn,:,:);
% %      aLDS(do(n), 1:NN,1:N) = LDS(do(n),:,:)-CLDS(mn,:,:);
% %      aLUS(do(n), 1:NN,1:N) = LUS(do(n),:,:)-CLUS(mn,:,:);
% %     end
% % 
% % 
% % end
% 
% aSWABS = aDT-aUT-aDS+aUS;
% %plot anomalies in global means
% x=sind(lat); 
% norm=trapz(x, ones(size(x))); %normalization for global means
% for nt=1:NT;
%     GDT(nt) =trapz(x, mean(squeeze(aDT(nt,:,:))') )/norm;
%     GUT(nt) =trapz(x, mean(squeeze(aUT(nt,:,:))') )/norm;
%     GOLR(nt) =trapz(x, mean(squeeze(aOLR(nt,:,:))') )/norm;
% 
%     GUS(nt) =trapz(x, mean(squeeze(aUS(nt,:,:))') )/norm;
%     GDS(nt) =trapz(x, mean(squeeze(aDS(nt,:,:))') )/norm;
%     GLDS(nt) =trapz(x, mean(squeeze(aLDS(nt,:,:))') )/norm;
%     GLUS(nt) = trapz(x, mean(squeeze(aLUS(nt,:,:))') )/norm;
%    
%     GSWABS(nt) = trapz(x, mean(squeeze(aSWABS(nt,:,:))') )/norm;
% 
% 
% end
% 
% 
% time=year+(month-0.5)/12;
% LWi=0.5; LW=2; LWA=3;
% FSz=28;
% FN= 'Arial';
% 
% cut_per=12; %cutoff period for low pass filter
% fs=1;
% 
% [low_GASR]=low_filter(GDT-GUT,fs, cut_per); 
% [low_GOLR]=low_filter(GOLR,fs, cut_per); 
% 
% [low_GDS]=low_filter(GDS,fs, cut_per); 
% [low_GUS]=low_filter(GUS,fs, cut_per); 
% [low_GLDS]=low_filter(GLDS,fs, cut_per); 
% [low_GLUS]=low_filter(GLUS,fs, cut_per); 
% 
% [low_GSWABS]=low_filter(GSWABS,fs, cut_per); 
% 
% 
% green=[0 .5 0];
% figure(1); clf;
% hold on; plot(time, GDT-GUT, 'r-', 'LineWidth', LWi)
% hold on; plot(time, low_GASR, 'r-', 'LineWidth', LW)
% 
% hold on; plot(time, GOLR, 'g-', 'Color', green,'LineWidth', LWi)
% hold on; plot(time, low_GOLR, 'g-', 'Color', green, 'LineWidth', LW)
% 
% [a_asr b_asr ci_asr R_asr] = lreg3(time, GDT-GUT);
% [a_olr b_olr ci_olr R_olr] = lreg3(time, GOLR);
% 
% hold on; plot(time, a_asr+b_asr*time, 'r--', 'LineWidth', LWA)
% hold on; plot(time, a_olr+b_olr*time, '--','Color', green, 'LineWidth', LWA)
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA);
% grid on; box on;
% xlabel('Year', 'FontSize', FSz, 'FontName', FN)
% ylabel('Global mean radiative Anomaly (W m^-^2)', 'FontSize', FSz, 'FontName', FN)
% title('Time Series of Global mean TOA radiative anomaly', 'FontSize', FSz, 'FontName', FN)
% 
% 
% figure(2); clf;
% 
% hold on; plot(time, low_GLUS, 'r-', 'LineWidth', LW)
% hold on; plot(time, low_GLDS, 'g-', 'Color', green,'LineWidth', LW)
% legend('LW UP', 'LW DOWN')
% hold on; plot(time, GLUS, 'r-', 'LineWidth', LWi)
% hold on; plot(time, GLDS, 'g-', 'Color', green,'LineWidth', LWi)
% 
% [a_us b_us ci_us R_us] = lreg3(time, GUS);
% [a_ds b_ds ci_ds R_ds] = lreg3(time, GDS);
% [a_lds b_lds ci_lds R_lds] = lreg3(time, GLDS);
% [a_lus b_lus ci_lus R_lus] = lreg3(time, GLUS);
% [a_swabs b_swabs ci_swabs R_swabs] = lreg3(time, GSWABS);
% 
% hold on; plot(time, a_lus+b_lus*time, 'r--', 'LineWidth', LWA)
% hold on; plot(time, a_lds+b_lds*time, '--','Color', green, 'LineWidth', LWA)
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA);
% grid on; box on;
% xlabel('Year', 'FontSize', FSz, 'FontName', FN)
% ylabel('Global mean surface radiative Anomaly (W m^-^2)', 'FontSize', FSz, 'FontName', FN)
% title('Time Series of Global mean surface LW radiative anomaly', 'FontSize', FSz, 'FontName', FN)
% 
% figure(3); clf;
% 
% hold on; plot(time, low_GUS, 'r-', 'LineWidth', LW)
% hold on; plot(time, low_GDS, 'g-', 'Color', green,'LineWidth', LW)
% legend('SW UP', 'SW DOWN')
% hold on; plot(time, GUS, 'r-', 'LineWidth', LWi)
% hold on; plot(time, GDS, 'g-', 'Color', green,'LineWidth', LWi)
% 
% 
% hold on; plot(time, a_us+b_us*time, 'r--', 'LineWidth', LWA)
% hold on; plot(time, a_ds+b_ds*time, '--','Color', green, 'LineWidth', LWA)
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA);
% grid on; box on;
% xlabel('Year', 'FontSize', FSz, 'FontName', FN)
% ylabel('Global mean surface radiative Anomaly (W m^-^2)', 'FontSize', FSz, 'FontName', FN)
% title('Time Series of Global mean surface SW radiative anomaly', 'FontSize', FSz, 'FontName', FN)
% 
% 
% figure(4); clf;
% 
% hold on; plot(time, low_GSWABS, 'r-', 'LineWidth', LW)
% %hold on; plot(time, low_GDS, 'g-', 'Color', green,'LineWidth', LW)
% legend('SW UP', 'SW DOWN')
% hold on; plot(time, GSWABS, 'r-', 'LineWidth', LWi)
% %hold on; plot(time, GDS, 'g-', 'Color', green,'LineWidth', LWi)
% 
% 
% hold on; plot(time, a_swabs+b_swabs*time, 'r--', 'LineWidth', LWA)
% %hold on; plot(time, a_ds+b_ds*time, '--','Color', green, 'LineWidth', LWA)
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA);
% grid on; box on;
% xlabel('Year', 'FontSize', FSz, 'FontName', FN)
% ylabel('Global mean surface radiative Anomaly (W m^-^2)', 'FontSize', FSz, 'FontName', FN)
% title('Time Series of Global mean surface SW radiative anomaly', 'FontSize', FSz, 'FontName', FN)


%zonal mean anomalies and trends
ZASR =  squeeze(mean(aDT-aUT,3));
ZOLR =  squeeze(mean(aOLR,3));
ZSWABS =  squeeze(mean(aSWABS,3));

for nn=1:NN;
    fl=squeeze(ZASR(:,nn));
    ZASR_low(1:NT, nn) = low_filter(fl, fs, cut_per);
    [a b ci R] = lreg3(time, fl);
    b_ASR(nn) = b;

    fl=squeeze(ZOLR(:,nn));
    ZOLR_low(1:NT, nn) = low_filter(fl, fs, cut_per);
    [a b ci R] = lreg3(time, fl);
    b_OLR(nn) = b;

    fl=squeeze(ZSWABS(:,nn));
    ZSWABS_low(1:NT, nn) = low_filter(fl, fs, cut_per);
    [a b ci R] = lreg3(time, fl);
    b_SWABS(nn) = b;

    fl=squeeze(ZUS(:,nn));
    ZUS_low(1:NT, nn) = low_filter(fl, fs, cut_per);
    [a b ci R] = lreg3(time, fl);
    b_SWABS(nn) = b;
end

[TM LAT] = meshgrid(time, lat);
colortrans;
conts = [-2*cmx:2:-2 2:2:2*cmx];
cmx=10;
figure(11); clf;
pcolor(TM, LAT, ZASR_low'); colorbar; 
hold on; contour(TM, LAT, ZOLR_low', conts, 'LineWidth', 3); colorbar; 
colormap(cmap); caxis([cmx*[-1 1]]); colorbar('LineWidth',3);

figure(12);clf;
plot(b_ASR, lat, 'r-', 'LineWidth', LW);
hold on; plot(b_SWABS, lat, 'r--', 'LineWidth', LW);
hold on; plot(b_OLR, lat, 'g-','Color', green, 'LineWidth', LW)


figure(13); clf;
pcolor(TM, LAT, ZSWABS_low'); colorbar; 
%hold on; contour(TM, LAT, ZOLR_low', conts, 'LineWidth', 3); colorbar; 
colormap(cmap); caxis([cmx*[-1 1]]); colorbar('LineWidth',3);