%compare the seasonal heating from surface fluxes and direct solar
%absorption
clear all; close all;

sty=['b- '; 'b--'; 'c- ';'c--';'g- ';'g--';'r- ';'r--';'k- ';'k--';'m- ';'m--'];

%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
cmap=[0 0 1;0 0 1; cmap];
MS=18;

CN=length(cmap);
figure(1); clf;
FSz=22; LW=2.8; FN='Arial'; %figure plottinf options
LW2=1.5; LWA=2.4; %for the phasing info
LWavg=1.5;
depths=[2.4 6 12 24 50]; NR=length(depths);
dol1=23; %levels to plot up
dol2=10;

FSz=22; 
FN='Arial';

xticl= [0 30 60 90];
xtic=sind(xticl);
xticl2= [-90 -60 -30 0 30 60 90];
xtic2=sind(xticl2);
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

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

     SHF = squeeze(mean(shiftdim(FS -US + DS,2))); %the surface heat flux upward from the ocean (not counting the solar that passes through the atmosphere

     NN=length(lat); N=length(lon); %dimensionality
     
     %define the absorbed radiation
     SWABS = squeeze(mean(shiftdim(DT - UT -DS + US,2))); %SW radiation absorbed within the atmosphere
     ZDT=squeeze(mean(shiftdim(DT,2)));
     ZUT=squeeze(mean(shiftdim(UT,2)));
     ALB=ZUT./ZDT;
        
     %plot the annual mean (non-solar weighted) albedo
     X=sind(lat);
     xmn=0; xmx=1;
     figure(1);hold on;
     plot(X,nanmean(ALB), 'Color', cmap2(cn,:), 'LineWidth', LW);
     hold on;  plot(X,mean(ZUT)./mean(ZDT), '--','Color', cmap2(cn,:), 'LineWidth', LW);
        %average poleward of a certaon latitude
        x=sind(lat); do=find(lat>=-30 & lat<=30); %sine of latitude (for area weighting over the polar cap
        
        xdo=x(do);
        norm=trapz(xdo,ones(size(xdo))); %normaliztion for averages over this region
     Xtic(nr) = cn;
     Xticl(nr)=depths(nr);
     
     ALB_bar(nr)= trapz(xdo, nanmean(ALB(:,do)))/norm; %spatial average of the temporal mean albedo 
     hold on;  plot(X,ALB_bar(nr)*ones(size(X)), '--','Color', cmap2(cn,:), 'LineWidth', LWavg);
     
     GALB(nr)= (trapz(xdo, squeeze(mean(ZUT(:,do))))/norm)/(trapz(xdo, squeeze(mean(ZDT(:,do))))/norm);
     GDT(nr)= (trapz(xdo, squeeze(mean(ZDT(:,do))))/norm);
     GUT(nr)= (trapz(xdo, squeeze(mean(ZUT(:,do))))/norm);
     hold on; plot(xmn, GALB(nr), '+', 'MarkerSize', MS, 'LineWidth', 3, 'Color', cmap2(cn,:));
     hold on; plot(xmx, GALB(nr), '+', 'MarkerSize', MS, 'LineWidth', 3, 'Color', cmap2(cn,:));
     

     S_bar_ALB_bar(nr)=trapz(xdo, squeeze(mean(ZDT(:,do))).*(1-squeeze(nanmean(ALB(:,do)))))/norm; %time mean solar and albedo
     
     MN1=4;
     MN2=10;
     figure(2);hold on; 
     hold on; plot(X,squeeze(ALB(MN1,:)), '-','Color', cmap2(cn,:), 'LineWidth', LW);
     hold on; plot(X,squeeze(ALB(MN2,:)), '--','Color', cmap2(cn,:), 'LineWidth', LW);

     %plot times eries of the albedo at a given altitde
     
     ILAT=find2(lat,20);
     figure(3); hold on; 
     T=1:12;
     plot(T,squeeze(ALB(:,ILAT)), '-', 'Color', cmap2(cn,:), 'LineWidth', LW);
     
     %%%%%%%%%%%%%%%% the latitudinal profile of albedo %%%%
     MNS=1:2:12;
     for mn=1:length(MNS);
         alb_sc=1.5;
             INS=squeeze(ZDT(mn,:));
         INS=.4*(INS)/mean(INS);
         figure(10+nr);
         hold on; plot(X, mn*ones(1,NN)+ alb_sc*ALB(MNS(mn),:), 'b', 'LineWidth', LW);
          hold on; plot(X, mn*ones(1,NN)+INS, 'r--', 'LineWidth', LW);
         hold on; plot(X, mn*ones(1,NN), '-', 'Color', [.4 .4 .5], 'LineWidth', LWA);
        
     
         
         
         figure(10+nr);hold on;
     end
     legend('Planetary Albedo', 'Insolation', 'Zero Axis');
     
     title( [num2str(depths(nr)) ' Slab depth seasonal cycle of planetary albedo'], 'FontSize', FSz, 'FontName',FN);
    xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
    ylabel('Month/ Planetary Albedo', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 3, 'XTick', xtic2, 'XTickLabel', xticl2, 'YTick', 1:length(MNS), 'YTickLabel', MNS)
axis([-1 1 0 length(MNS)+1]);
hold on; plot([0 0], [0 length(MNS)+1], 'k--', 'LineWidth', 2);
         
     %%%%%%%%%%%%%%%% the latitudinal profile of albedo %%%%
     MNS=1:1:12;
     for mn=1:length(MNS);
         alb_sc=1.5;
             INS=squeeze(ZDT(mn,:));
         INS=.4*(INS)/mean(INS);
         figure(20+nr);
         hold on; plot(X, alb_sc*ALB(MNS(mn),:), char(sty(MNS(mn),:)), 'LineWidth', LW);
          
         %hold on; plot(X, mn*ones(1,NN)+INS, 'r--', 'LineWidth', LW);
         %hold on; plot(X, mn*ones(1,NN), '-', 'Color', [.4 .4 .5], 'LineWidth', LWA);
        
     
         
         
        
     end
     %legend('Planetary Albedo', 'Insolation', 'Zero Axis');
     
     title( [num2str(depths(nr)) ' Slab depth seasonal cycle of planetary albedo'], 'FontSize', FSz, 'FontName',FN);
    xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
    ylabel('Planetary Albedo', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 3, 'XTick', xtic2, 'XTickLabel', xticl2, 'YTick', 1:length(MNS), 'YTickLabel', MNS)
ymn=.2; ymx=0.8;
axis([-1 1 ymn ymx]);
hold on; plot([0 0], [ymn ymx], 'k--', 'LineWidth', 2) ;
legend('JAN','FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL','AUG', 'SEP','OCT', 'NOV','DEC')

     
     
end




     
cut_lat=0;
figure(1);hold on;
colormap(cmap2);colorbar;
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
ylabel('Time Mean Planetary Albedo', 'FontSize', FSz, 'FontName', FN);
title('Annual mean plantery albedo (not solar weighted)', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 3, 'XTick', xtic, 'XTickLabel', xticl)
H=colorbar('location', 'SouthOutside', 'XTick',Xtic, 'XTickLabel', Xticl, 'LineWidth', 3, 'FontSize', FSz, 'FontName', FN); 
axis([xmn xmx  0 .55]); box on;  

figure(2);hold on;
colormap(cmap2);colorbar;
xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
ylabel('Time Mean Planetary Albedo', 'FontSize', FSz, 'FontName', FN);
title('Annual mean plantery albedo (not solar weighted)', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 3, 'XTick', xtic, 'XTickLabel', xticl)
H=colorbar('location', 'SouthOutside', 'XTick',Xtic, 'XTickLabel', Xticl, 'LineWidth', 3, 'FontSize', FSz, 'FontName', FN); 
axis([-1 xmx  0 .55]); box on; 

%% TIME SERIES %%%%%%%%%%%%%
INS=squeeze(DT(:,ILAT));
INS=.3+.1*(INS-mean(INS))/std(INS);
figure(3);hold on; plot(T, INS, 'k--', 'LineWidth',3);
hold on; colormap(cmap2);colorbar;
xlabel('Month', 'FontSize', FSz, 'FontName', FN);
title(['Planetary Albedo albeod at ' num2str(lat(ILAT)) '^o latitude'], 'FontSize', FSz, 'FontName', FN);
ylabel({'Planetary Albedo and';'insolation in nonsense units'}, 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 3)
H=colorbar('location', 'SouthOutside', 'XTick',Xtic, 'XTickLabel', Xticl, 'LineWidth', 3, 'FontSize', FSz, 'FontName', FN); 


