%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;

     
     cut_lat=30;
      solstice=174; %where the extratropical DT phase is
     FSz=20; FN='Arial';
     LWA=2.5;
     ytic=-900:200:-100;
     yticl=-ytic;

     %%% 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);

      xmn=-90; xmx = 90; ymn=-1000; ymx=0;
      figure(1);clf;
      MN=20; MX=60;
      REF=NaN*ones(NL,NN);
      REF(1:NL, find(tlat>=30))=solstice;
      solstice_sh=solstice-365/2;
      REF(1:NL, find(tlat<=-30))=solstice_sh;
      FIELD=squeeze(OBS_PHASE)-REF;
      FIELD(find(FIELD >= 365)) =  FIELD(find(FIELD >= 365)) -365; 
      pcolor(LAT, -LEV,FIELD);
      colormap('jet'); caxis([MN MX]); colorbar;
      shading interp;
      axis([xmn xmx ymn ymx])
      set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'YTick', ytic, 'YTickLabel', yticl)
      xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
      ylabel('Pressure Level (hPa)', 'FontSize', FSz, 'FontName', FN);
      title([ 'OBSERVED phase'], 'FontSize', FSz, 'FontName', FN)
      box on;
      
      %now repeat with the PI models
      clear REF;
      load('C:\Users\thedhoe\Desktop\aqua_depth\figures\phase_obs_CMIP3\PI_T_common','CRNM','C_T','clon', 'clat','cplev'); %compress the zonal dimension -- is already averaged
          
      NAMES={'GFDL - AM2.0';
    'BCCR';
    'CAM3';
    'cccma';
    'CNRM';
    'CSIRO';
    'ECHAM';
    'ECHO';
    'HADCM3';
    'INMCM3';
    'IPSL';
    'MRI';}
      
      [NM, NT, NL, NN, N] = size(C_T);
      [LAT, LEV] = meshgrid(clat, cplev);
      for nm=1:NM;
          ZT = squeeze(nanmean(shiftdim( squeeze( C_T(nm,:,:,:,:)), 3)));
          ph=fund_phase_3d(ZT);
          REF=NaN*ones(NL,NN);
          REF(1:NL, find(clat>=30))=solstice;
          solstice_sh=solstice-365/2;
          REF(1:NL, find(clat<=-30))=solstice_sh;
          FIELD=squeeze(ph)-REF;
          FIELD(find(FIELD >= 365)) =  FIELD(find(FIELD >= 365)) -365; 
          figure(10+nm); clf;
          pcolor(LAT, -LEV/100,FIELD);
          colormap('jet'); caxis([MN MX]); colorbar;
          shading interp;
          axis([xmn xmx ymn ymx])
          set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'YTick', ytic, 'YTickLabel', yticl)
          xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
          ylabel('Pressure Level (hPa)', 'FontSize', FSz, 'FontName', FN);
          title([ char(NAMES(nm,:)) ' phase'], 'FontSize', FSz, 'FontName', FN)
          box on;
      end
          
          
      % and the ensemble mean
          ZT = squeeze(nanmean(shiftdim( squeeze( nanmean(C_T(:,:,:,:,:))), 3)));
          ph=fund_phase_3d(ZT);
          REF=NaN*ones(NL,NN);
          REF(1:NL, find(clat>=30))=solstice;
          solstice_sh=solstice-365/2;
          REF(1:NL, find(clat<=-30))=solstice_sh;
          FIELD=squeeze(ph)-REF;
          FIELD(find(FIELD >= 365)) =  FIELD(find(FIELD >= 365)) -365; 
          figure(2); clf;
          pcolor(LAT, -LEV/100,FIELD);
          colormap('jet'); caxis([MN MX]); colorbar;
          shading interp;
          axis([xmn xmx ymn ymx])
          set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', LWA, 'YTick', ytic, 'YTickLabel', yticl)
          xlabel('Latitude', 'FontSize', FSz, 'FontName', FN);
          ylabel('Pressure Level (hPa)', 'FontSize', FSz, 'FontName', FN);
          title(['Ensemble average phase'], 'FontSize', FSz, 'FontName', FN)
          box on;
      
      

      
  
      