% %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=174; %where the extratropical DT phase is
%      
% 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];
% 
% 
 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']);
%      NL=length(Pfull); NN=length(lat); N=length(lon);
%      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
%      
%      for nl=1:NL;
%          for nn=1:NN;
%             phase_cross(nr,nl,nn)=fund_phase(squeeze(ZT(:,nl,nn)));
%          end
%      end
%      
      xmn=30; xmx=90;
      ymn=-1000; ymx=-300;
      MN=45; MX=115;
%      
     ytic=ymn:100:ymx;
      yticl=-ytic;
%      [LAT LEV]=meshgrid(lat, Pfull);
%      figure(nr);
%      pcolor(LAT, -LEV, squeeze(phase_cross(nr,:,:))-solstice);
%      colormap('jet'); caxis([MN MX]); colorbar;
%      shading interp;
%      axis([xmn xmx ymn ymx])
%      set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTick', ytic, 'YTickLabel', yticl)
%      xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
%      ylabel('Pressure Level (hPa)', 'FontSize', FSz, 'FontName', FN);
%      title([num2str(depths(:, nr)) ' meters mixed layer'], 'FontSize', FSz, 'FontName', FN)
%      
%      
%      
% end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OBSERVATIONS %%%%%%%%%%%%%%%%%%%%%%%%%
%      
%      %%% add in the observations
%      
%       load C:\Users\thedhoe\Desktop\CMIP_seasonal_lap\obs\data\T.mat
%       tlat=lat;
%       tlon=lon;
%       [TLON, TLAT]=meshgrid(tlon, tlat);
%       [NT, NL, NN, N]=size(T);
%       
%       load('C:\Users\thedhoe\Desktop\aqua_depth\figures\phase_cross_sections\OBS_mean_PS', 'lat', 'lon','PS')
%       pslon=double(lon);
%       pslat=double(lat);
%       
%       [PSLON, PSLAT]=meshgrid(pslon, pslat);
%       
%       %interpolate the surface pressure to the temperature gbrid
%       IPS=interp2(PSLON, PSLAT,PS, TLON, TLAT);
%       
%       MASK=ones(NL,NN,N);
%       
%       
%       for nn=1:NN;
%           for n=1:N;
%               do=find(plev >= (squeeze(IPS(nn,n)))/100); %levels that are below the surface at the given grid point
%               MASK(do,nn,n)=NaN; %set these levels equal to NaN in the Mask
%           end
%       end
      
      %compose a zonal mean-- removing the gridpoints that are below the
      %surface
      for mn=1:12; %loop over month
        ZT(mn,1:NL, 1:NN)=squeeze(nanmean(shiftdim( squeeze(T(mn,:,:,:)).*MASK,2)));
      end
      
      %take the phase of the zonal mean temperature
      for nn=1:NN;
              for nl=1:NL;
                  
                   OBS_PHASE(nl,nn) = fund_phase(squeeze(ZT(:,nl,nn)));
                   OBS_AMP(nl,nn) = fund(squeeze(ZT(:,nl,nn)));
              end
      end
      
      
      [LAT, LEV]=meshgrid(tlat,plev);
      MN=0; MX=70;
      
      figure(19);clf;
      pcolor(LAT, -LEV,squeeze(OBS_PHASE)-solstice);
      colormap('jet'); caxis([MN MX]); colorbar;
      shading interp;
      axis([xmn xmx ymn ymx])
      set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTick', ytic, 'YTickLabel', yticl)
      xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
      ylabel('Pressure Level (hPa)', 'FontSize', FSz, 'FontName', FN);
      title([ 'OBSERVED phase -- NH '], 'FontSize', FSz, 'FontName', FN)
      
      xmn_sh=-90; xmx_sh=-30;
      solstice_sh=solstice-365/2;
      figure(20);clf;
      FIELD=squeeze(OBS_PHASE)-solstice_sh;
      FIELD(find(FIELD > 365/2)) = FIELD(find(FIELD > 365/2)) -365;
      FIELD(find(FIELD < 0)) = FIELD(find(FIELD < 0)) +365;
      pcolor(LAT, -LEV,FIELD);
      colormap('jet'); caxis([MN MX]); colorbar;
      shading interp;
      axis([xmn_sh xmx_sh ymn ymx])
      set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2, 'YTick', ytic, 'YTickLabel', yticl)
      xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
      ylabel('Pressure Level (hPa)', 'FontSize', FSz, 'FontName', FN);
      title([ 'OBSERVED phase -- NH '], 'FontSize', FSz, 'FontName', FN)
      
      