%look at the global mean temperature difference and the energy budget

clear all; close all;
%load the climatologiclay data
dr =    {'C:/research/aqua_gcm/slab_2.4m/';  %directory listing dfor each run
         'C:/research/aqua_gcm/slab_6m/';
         'C:/research/aqua_gcm/slab_12m/';
         'C:/research/aqua_gcm/slab_24m/';
         'C:/research/aqua_gcm/slab_50m/'};
     
NR=5;
colormap('jet');cmap=colormap; %save the RGB colormap vector
CN=length(cmap);
figure(1); clf;
FSz=22; LW=2; FN='Arial'; %figure plottinf options
LW2=1.5; %for the phasing info
off=10; %opffset to chop off the end of the colobar
for n=1:3
cmap2(1:CN-off,n)=interp1((0:(CN-1-off))/(CN-1-off), cmap(1:end-off,n), ((0:(CN-1-off))/(CN-1-off)).^.4); %elongate the shallow end of the color scale to add contrast there
end

cross=30;
depths=[2.4 6 12 24 50];
for nr=1:5;  %loop over run (depths)
     cn=round((CN-off)*depths(nr)/depths(NR)); %index of the given run in the colomap matrix (could make it scale as depth if you want
     load([char(dr(nr)) 'energy_climatology.mat']);
     load([char(dr(nr)) 'TV_trop.mat']);
     load([char(dr(nr)) 'T.mat']);
     
     tdo=find(lat >=-cross & lat <=cross); 
     edo=find(lat >=cross);
     
     ATS=mean(squeeze(mean(squeeze(T(:,end,:,:))))'); %zonal annual mean surface temperature
     AALB=mean(squeeze(mean(squeeze(UT(:,:,:))))')./mean(squeeze(mean(squeeze(DT(:,:,:))))');
     TV=TV_trop;
     FS = FS -US + DS; %the surface heat flux upward from the ocean (not counting the solar that passes through the atmosphere

     NN=length(lat); N=length(lon); %dimensionality
     NL=length(Pfull);
     %define the absorbed radiation
     absorb = DT - UT -DS + US; %SW radiation absorbed within the atmosphere
     Zabsorb=squeeze(mean(shiftdim(absorb,2))); %zonal mean absorption
     ZDT=squeeze(mean(shiftdim(DT,2)));
     ZUT=squeeze(mean(shiftdim(UT,2)));
     ZFS=squeeze(mean(shiftdim(FS,2)));
     X=sind(lat);
     NORM=trapz(X,ones(size(X))); %normalization for global mean
     GTS(nr)=trapz(X, ATS)/NORM;
     GDT(nr)=trapz(X,squeeze(mean(ZDT)))/NORM;
     GUT(nr)=trapz(X,squeeze(mean(ZUT)))/NORM;
     
     %partition the upwelling solar at TOA into tropical and extratropical
     %pieces
     
     TUT(nr)=trapz(X(tdo), squeeze(mean(ZUT(:,tdo))))/NORM;
     EUT(nr)=2*trapz(X(edo), squeeze(mean(ZUT(:,edo))))/NORM;
     xx=sind(lat);
     figure(19); hold on; 
     plot(xx,ATS, 'Color', cmap2(cn,:), 'LineWidth', LW)
     %hold on; plot(-lat, ATS,'Color', cmap2(cn,:),'LineStyle', '--', 'LineWidth', LW) 
  
     figure(20); hold on; 
     plot(xx,AALB, 'Color', cmap2(cn,:), 'LineWidth', LW)

     
     ANN_absorb(nr,1:NN) = mean(squeeze(mean(absorb))');
     ANN_FS(nr, 1:NN) = mean(squeeze(mean(FS))');
     
     Xtic(nr) = cn;
     Xticl(nr)=depths(nr);
end

xticl=-90:30:90;
xtic=sind(xticl);
cut_lat=0;
ymn=260;
ymx=310;
xmn=sind(cut_lat);
xmx=sind(90);
figure(19); hold on;
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
ylabel('Annual Mean Surface Temperature', 'FontSize', FSz, 'FontName', FN)
H=colorbar('location', 'SouthOutside', 'XTick',Xtic, 'XTickLabel', Xticl, 'LineWidth', 3, 'FontSize', FSz, 'FontName', FN); 
set(gca,'FontSize', FSz, 'FontName', FN, 'LineWidth', 2.5, 'XTick', xtic, 'XTickLabel', xticl)
box on;
axis([xmn xmx ymn ymx]);
     
ymn=0.1;
ymx=0.5;
figure(20); hold on;
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
ylabel('Annual Mean Planetary Albedo', 'FontSize', FSz, 'FontName', FN)
H=colorbar('location', 'SouthOutside', 'XTick',Xtic, 'XTickLabel', Xticl, 'LineWidth', 3, 'FontSize', FSz, 'FontName', FN); 
set(gca,'FontSize', FSz, 'FontName', FN, 'LineWidth', 2.5, 'XTick', xtic, 'XTickLabel', xticl)
box on;
axis([xmn xmx ymn ymx]);
     

%plot as temperature anomaly relative to
LWA=2.5;
REF=5; %run to reference to
BW=0.5;
BW2=0.3;
LWB=2;
T=GTS-GTS(REF);
NET=-(GUT-GUT(REF));
sc=0.5; %scaling between radiation and temperature
pos=1:2:9; off=0.5;
xtic=pos;
xticl={'2.4m'; '6m'; '12m';'24m';'50m'}; 
figure(23); clf;
hold on; bar(pos,T, 'b', 'BarWidth', BW, 'LineWidth', LWB)
hold on; bar(pos+off,sc*NET, 'r', 'BarWidth', BW2, 'LineWidth', LWB)

ymn=-10; ymx=5;
ymn2=-20;
ymx2=10;
xmn=0; xmx=10;
yticl2=ymn2:10:ymx2;
%ytic2=(yticl2/sc)
ax(1)=gca;
set(ax(1),'Position',[0.12 0.12 0.75 0.70])
set(ax(1),'XColor','k','YColor','k', 'LineWidth', LWA);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'XTick', xtic, 'XTickLabel', xticl)

axis([xmn xmx ymn ymx]);
ax(2)=gca;
%set(ax(2),'Position',[0.12 0.12 0.75 0.70]);
ax(2)=axes('Position',get(ax(1),'Position'),...
   'XAxisLocation','top',...
   'YAxisLocation','right',...
   'Color','none',...
   'XColor','k','YColor','k');
set(ax(2),'XColor','k','YColor','r', 'LineWidth', LWA);
axis([xmn xmx ymn2 ymx2]);
ylabel('Radiation anomaly (relative to 50 m run in W m^-^2)', 'FontSize', FSz, 'FontName', FN)
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'XTick', xtic, 'XTickLabel', xticl);
%data=SF*[cNET_tot cSW_tot -cLW_tot; cNET_trop cSW_trop -cLW_trop;cNET_sub cSW_sub -cLW_sub; cNET_ex cSW_ex -cLW_ex];

% figure(7);clf;
% space=.7; sep=.12;
% %h=bar([-sep sep  space-sep space+sep 2*space-sep 2*space+sep  3*space-sep 3*space+sep],data, 'grouped','BarWidth', 1);
% h=bar([0  space 2*space 3*space],data, 'grouped','BarWidth', 1);
% xtic=[0 space 2*space 3*space];
% xticl={'Total'; 'Tropics';'Subtropics';'Extratropics'};
% ylabel('Contribution to MHT_E_Q (PW)', 'FontSize', FSz, 'FontName', FN)
% set(gca, 'FontSize', FSz, 'FontName', FN, 'XTick', xtic, 'XTickLabel', xticl)
% title('CERES Hemispheric asymetry in TOA radiation by region', 'FontSize', FSz, 'FontName', FN)
% colormap([0 0 0; 1 0 0; 0 1 0]);
% box on;
% legend('Net TOA radiation', 'SW', 'LW')