%look at the amplitude and phase of the tempearure variations 
%compare the seasonal heating from surface fluxes and direct solar
%absorption
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/'};
     
     cut_lat=30;
     solstice=monthday(6)+21;
     
NR=5;
colormap('jet');cmap=colormap; %save the RGB colormap vector
cmap=[0 0 1;0 0 1; cmap];

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

FSz=22; 
FN='Arial';


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
     %cn=round(CN*(nr)/(NR));
     load([char(dr(nr)) 'T.mat']);
     ZT=squeeze(mean(shiftdim(T,3))); %zonal mean Temperature (month, pressure_levl, lat)
     mlev=[0; mean([Pfull(1:end-1,1) Pfull(2:end,1)]')'; 1000];
     weight=diff(mlev); weight(find(mlev<400))=0; %mass weighted vector for the troposphere only
     norm=weight/sum(weight); %normalized weighting index
     
     
     do=find(lat>=cut_lat);
     X=lat; NN=length(lat);
     NL=length(Pfull);
     for mn=1:12;
         VZT(mn,1:NN)=norm'*squeeze(ZT(mn,:,:));
     end
     
     AVZT = fund(VZT');
     AZT=fund_3d(ZT);
     
     for nn=1:NN;
         PHASE(nr,nn) = fund_phase(squeeze(VZT(:,nn)));
         AMP(nr,nn) = fund(squeeze(VZT(:,nn)));
         
         for nl=1:NL;
          phase(nl,nn) = fund_phase(squeeze(ZT(:,nl,nn)));
          
         end
     end
     
     
     figure(1); hold on;
     plot(X(do),AZT(dol1,do),'-','Color', cmap2(cn,1:3), 'LineWidth', LW);
     hold on; plot(X(do),AZT(dol2,do),'--','Color', cmap2(cn,1:3), 'LineWidth', LW)
     
     figure(2); hold on;
     plot(X(do),phase(dol1,do)-solstice,'-','Color', cmap2(cn,1:3), 'LineWidth', LW);
     hold on; plot(X(do),phase(dol2,do)-solstice,'--','Color', cmap2(cn,1:3), 'LineWidth', LW)
     
     figure(3); hold on;
     plot(X(do),PHASE(nr,do)-solstice,'-','Color', cmap2(cn,1:3), 'LineWidth', LW);
    
     
     Xtic(nr) = cn;
     Xticl(nr)=depths(nr);
     
end
colormap(cmap2);
figure(1); hold on;
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 3)
 xlabel('Latitude', 'FontSize', FSz', 'FontName', FN)
 ylabel('Seasonal Amplitude of temperature (K)','FontSize', FSz, 'FontName', FN) 
 H=colorbar('location', 'SouthOutside', 'XTick',Xtic, 'XTickLabel', Xticl, 'LineWidth', 3, 'FontSize', FSz, 'FontName', FN); 
 axis([cut_lat 90 0 45])
 box on;    

figure(2); hold on;

     xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
     ylabel('Seasonal phase lag of temperature (days relative to solstice)', '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)
     box on;
     axis([cut_lat 90 50 100])
     
     figure(3); hold on;

     xlabel('Latitude', 'FontSize', FSz, 'FontName', FN)
     ylabel('Seasonal phase lag of temperature (days relative to solstice)', '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)
     box on;
     axis([cut_lat 90 50 100])
     
     %%% add in the observations
     
      load C:\Users\thedhoe\Desktop\CMIP_seasonal_lap\obs\data\T.mat
 elat=lat; elon=lon; %save the ERA lat and lon
load C:\Users\thedhoe\Desktop\CMIP_seasonal_lap\obs\data\MASK.mat
load('C:\Users\thedhoe\Desktop\CMIP_seasonal_lap\obs\data\cross_mask.mat', 'cross_mask', 'clon', 'clat', 'cplev') %the topography mask
[CLAT, CLEV]=meshgrid( [-90 clat 90], [cplev/100 0]); %extend the cross section topography mask
cf=[cross_mask(:,1) cross_mask cross_mask(:,end)]; %extend data latitude
cf=[cf; cf(end,:)]; %extend data longtidue

[NT, NL, blah, blah2]=size(T);

%plot the seasonal amplitudes of temperature
SA_T=fund_3d(squeeze(mean(shiftdim(T,3)))); %seasonal amplitude of zonal mean temperature (cross section)

do1_obs=length(plev);
do2_obs=find2(plev,600);
do_NH=find(elat>=cut_lat);
do_SH=find(elat<=-cut_lat);



ZT=squeeze(mean(shiftdim(T,3))); %zonal mean
AZT=fund_3d(ZT); %seasonal amplitude


for nn=1:length(elat);
    for nl=1:length(plev);
        obs_phase(nl,nn)=fund_phase(squeeze(ZT(:,nl,nn)));
    end
end

figure(1); hold on;
     hold on; plot(elat(do_NH),AZT(do1_obs,do_NH),'k-', 'LineWidth', LW);
     hold on; plot(elat(do_NH),AZT(do2_obs,do_NH),'k--', 'LineWidth', LW);
     hold on; plot(-elat(do_SH),AZT(do1_obs,do_SH),'m-',  'LineWidth', LW);
     hold on; plot(-elat(do_SH),AZT(do2_obs,do_SH),'m--', 'LineWidth', LW);
     
     figure(2); hold on;
     hold on; plot(elat(do_NH),obs_phase(do1_obs,do_NH)-solstice,'k-', 'LineWidth', LW);
     hold on; plot(elat(do_NH),obs_phase(do2_obs,do_NH)-solstice,'k--', 'LineWidth', LW);
%      hold on; plot(-elat(do_SH),obs_phase(do1_obs,do_SH)-solstice,'m-',  'LineWidth', LW);
%      hold on; plot(-elat(do_SH),obs_phase(do2_obs,do_SH)-solstice,'m--', 'LineWidth', LW);
 