% make composite zonal mean temperature maps of:
%     A: annual mean anomaly from global mean temp at each level
%     B: Seasonal cycle (amplitude and phase)

clear all; close all;
%load the 4-d maps od temperature from each aqyaplanet run
dr =    {'/home/disk/viento/aaron/research/aqua_gcm/slab_2.4m/';
         '/home/disk/viento/aaron/research/aqua_gcm/slab_6m/';
         '/home/disk/viento/aaron/research/aqua_gcm/slab_12m/';
         '/home/disk/viento/aaron/research/aqua_gcm/slab_24m/';
         '/home/disk/viento/aaron/research/aqua_gcm/slab_50m/'};

     FSz=21; FN='Arial';
     
     CMAX=22;
     depths = [2.4 6 12 24 50];
for nr=1:5;

    %load the 4d temps
    load([char(dr(nr)) 'T.mat']);
    [NT, NL, NN, N] = size(T); %dimensionality
    x=sin(lat*pi/180); %sine of latitude
    norm=trapz(x, ones(size(x))); %normaliztion for global integrals (should equal 2)
    ZAT = squeeze(mean(shiftdim(squeeze(mean(T)),2))); %zonal annual mean temp (level,latitude)
    ZT = squeeze(mean(shiftdim(T,3))); %zonal mean temperature (month,level, latitude)
    
    for nl=1:NL; %loop over vertical level
        T_bar(nl) = trapz(x,ZAT(nl,:))/norm; %global annual mean at each level
        TAp(nl,1:NN) = ZAT(nl,:) - T_bar(nl); %annual mean spatial anomaly from global mean by level 
    
        [A, ph] = fund2(squeeze(ZT(:,nl,:))');
        SEAS(nr,nl,1:NN) = A;
        SEAS_phase (nr,nl,1:NN) = ph;
        
    end
    
    
    dol=find(Pfull >= 250); %tropopause cutoff for plotting;
    
    ytic=-1000:100:-200;
    yticl=char(-ytic);
    
    [LAT, LEV] = meshgrid(lat,Pfull); %mesh for zonal mean by level
    color_trans; %red and blue colorbar
    colormap(cmap);
    
    ANN(nr,1:NL,1:NN) = TAp; %save the annual mean anomaly field  
    
end
    
    %cmax=0.9*max(abs(g2c(Tp))); %maximum value of the temp anomaly field

figure(1); clf;
%subplot(3,2,1)
pcolor(LAT(dol,:), -LEV(dol,:), squeeze(mean(ANN(:,dol,:)))); 
set(gca, 'FontSize', FSz,'FontName', FN); %'YTick', ytic, 'YTickLabel', yticl);
xlabel('Latitude', 'FontSize',FSz,'FontName', FN);
ylabel('Pressure Level (hPA)','FontSize', FSz,'FontName', FN);
title('Ensemble Mean Spatial Anomaly', 'FontSize', FSz,'FontName', FN);
caxis([-CMAX CMAX]); colorbar;
set(gca, 'FontSize', FSz,'FontName', FN); %'YTick', ytic, 'YTickLabel', yticl);

CMAX2=10;

for nr=1:4:5   
    figure(nr+1); hold on;
    %subplot(3,2,nr+1)
    pcolor(LAT(dol,:), -LEV(dol,:), squeeze(ANN(nr,dol,:))-squeeze(mean(ANN(:,dol,:)))); 
    caxis([-CMAX2 CMAX2]); %center the colorbar scale
    colormap(cmap); colorbar;
    
    set(gca, 'FontSize', FSz,'FontName', FN); %'YTick', ytic, 'YTickLabel', yticl);
    xlabel('Latitude', 'FontSize',FSz,'FontName', FN);
    ylabel('Pressure Level (hPA)','FontSize', FSz,'FontName', FN);
    title([num2str(depths(nr)) ' Meter Anomaly from Ensemble mean'], 'FontSize', FSz, 'FontName', FN)
end

for nr=1:5;
        
    figure(nr+10); hold on;
    %subplot(3,2,nr+1)
    pcolor(LAT(dol,:), -LEV(dol,:), squeeze(SEAS(nr,dol,:))-0*squeeze(mean(SEAS(:,dol,:)))); 
    %caxis([-CMAX2 CMAX2]); %center the colorbar scale
    colormap('jet'); colorbar;
    
    set(gca, 'FontSize', FSz,'FontName', FN); %'YTick', ytic, 'YTickLabel', yticl);
    xlabel('Latitude', 'FontSize',FSz,'FontName', FN);
    ylabel('Pressure Level (hPA)','FontSize', FSz,'FontName', FN);
    title([num2str(depths(nr)) ' Seasonal Amplitude Temperature '], 'FontSize', FSz, 'FontName', FN)
end
        
cd figures