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]));
 
UTcs(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'toa_sw_clr_t_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
UTcs((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'toa_sw_clr_t_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])); 

OLRcs(1:NT1, 1:NN,1:N) = double(permute(ncread(fil1, 'toa_lw_clr_t_mon', [1 1 1], [Inf Inf Inf]), [3 2 1]));
OLRcs((NT1+1):NT, 1:NN,1:N) = double(permute(ncread(fil2, 'toa_lw_clr_t_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])); 

save('C:\Users\Aaron Donohoe\Desktop\research\radiation\ceres_ebaf\4.1\CERES_TOA.mat', 'lat', 'lon', 'year', 'month', 'OLR', 'DT', 'UT', 'OLRcs', 'UTcs')

return
%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));
ZUS =  squeeze(mean(aUS,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_US(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, ZUS_low', conts, 'LineWidth', 3); colorbar; 
colormap(cmap); caxis([cmx*[-1 1]]); colorbar('LineWidth',3);