%compute AHT (especially at the equator) in the TRACE experiment
clear all; close all;
load('F:\trace\LH_dec.mat') %Latent heat flux
load('F:\trace\SH_dec.mat') %Sensible heat flux
load('F:\trace\FLNS_dec.mat') %net longwave flux at the surface
load('F:\trace\FSNS_dec.mat') %net shortwave flux at the surface
load('F:\trace\OLR_dec.mat') %OLR - top of model
load('F:\trace\FSNT_dec.mat') %net shortwave at TOA


 %make the surface heat flux -- to the atmosphere
SHF= -FSNS + LH + SH + FLNS;
TOA = FSNT-OLR;
ZTOA = squeeze(mean(shiftdim(TOA,2)));
ZSHF = squeeze(mean(shiftdim(SHF,2)));
ZOLR = squeeze(mean(shiftdim(OLR,2)));
ZFSNT = squeeze(mean(shiftdim(FSNT,2)));

NY=length(YEARS);
NN=length(lat);
N=length(lon)
x=sind(lat); %sin of latitude as a corrdinate for global integrlas (or avverages)
ae=6.36*10^6; %Earth radius
norm=trapz(x, ones(size(x))); %normalization for spatial weighted global averages

 
 mlat=[-90; (lat(2:end)+lat(1:end-1))/2;90]; mNN=length(mlat); %the heat transport will be calculated on this grid 
 dlat=diff(mlat);
 dA=(2*pi*ae*cosd(lat)).*(ae*dlat*pi/180);
 
 for ny=1:NY; %for ny=1:NY;
     ny;
 GTOA(ny)=trapz(x, squeeze(ZTOA(ny,:)))/norm;
 GSHF(ny)=trapz(x, squeeze(ZSHF(ny,:)))/norm;
 GASR(ny)=trapz(x, squeeze(ZFSNT(ny,:)))/norm;
 GLW(ny)=trapz(x, squeeze(ZOLR(ny,:)))/norm;
 
 AHT(ny,1)=0; OHT(ny,1)=0;
 
 LW_HT(ny,1)=0; SW_HT(ny,1)=0;
 for nn=2:mNN;
 AHT(ny,nn) = AHT(ny,nn-1) + dA(nn-1)*(ZTOA(ny,nn-1) - GTOA(ny) + ZSHF(ny,nn-1)-GSHF(ny));
 OHT(ny,nn) = OHT(ny,nn-1) + dA(nn-1)*(-ZSHF(ny,nn-1)+GSHF(ny));
 LW_HT(ny,nn) = LW_HT(ny,nn-1) + dA(nn-1)*(-ZOLR(ny,nn-1)+GLW(ny));
 SW_HT(ny,nn) = SW_HT(ny,nn-1) + dA(nn-1)*(ZFSNT(ny,nn-1)-GASR(ny));
 end
 end
 

 
  %[low_AHT]=low_filter(AHT_EQ,1,cut_per); %200 year low pass filter
 
  FSz=28; FN='Arial';
  PW=10^-15; LW=2;
  figure(1);clf;
  hold on; plot(YEARS, PW*max(AHT'), 'r', 'LineWidth', LW);
  hold on; plot(YEARS, PW*max(AHT'+OHT'), 'm', 'LineWidth', LW);
  hold on; plot(YEARS, -PW*min(AHT'), 'b', 'LineWidth', LW);
  hold on; plot(YEARS, -PW*min(AHT'+OHT'), 'c', 'LineWidth', LW);
  ymx=6.6; ymn=4.8;
  YD=-14400;
  hold on; plot(YD*[1 1], [ymn ymx], 'k--', 'LineWidth', 2)
  
  ymx=3; ymn=11;
  figure(3); clf;
  hold on; plot(YEARS, PW*max(SW_HT'), 'r', 'LineWidth', LW)
  hold on; plot(YEARS, PW*max(-LW_HT'), 'g', 'LineWidth', LW)
  hold on; plot(YEARS, PW*max(AHT'+OHT'), 'k', 'LineWidth', LW)
  hold on; plot(YEARS, PW*( -max(-LW_HT') + max(SW_HT')) , 'm', 'LineWidth', LW)
  %check to see if they are additivie
  hold on; plot(YEARS, PW*( -max(-LW_HT') + max(SW_HT')) , 'm', 'LineWidth', LW)
  title('NH total heat transport decomposed into ASR* and OLR*', 'FontSize', FSz, 'FontName', FN)
  set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 3); box on;
  xlabel('Years before present', 'FontSize', FSz, 'FontName', FN);
  ylabel('MHT_M_A_X , ASR*, OLR* (PW)', 'FontSize', FSz, 'FontName', FN);
  legend('AHT_M_A_X', 'ASR*', 'OLR*', 'ASR*-OLR*')
  hold on; plot(YD*[1 1], [ymn ymx], 'k--', 'LineWidth', 2)
  
  figure(2); clf;
  hold on; plot(YEARS, -PW*min(SW_HT'), 'r', 'LineWidth', LW)
  hold on; plot(YEARS, PW*max(LW_HT'), 'g', 'LineWidth', LW)
  hold on; plot(YEARS, -PW*min(AHT'+OHT'), 'k', 'LineWidth', LW)
  %check to see if they are additivie
  hold on; plot(YEARS, PW*( -max(LW_HT') - min(SW_HT')) , 'm', 'LineWidth', LW)
  title('SH total heat transport decomposed into ASR* and OLR*', 'FontSize', FSz, 'FontName', FN)
  set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 3); box on;
  xlabel('Years before present', 'FontSize', FSz, 'FontName', FN);
  ylabel('MHT_M_A_X , ASR*, OLR* (PW)', 'FontSize', FSz, 'FontName', FN);
  legend('AHT_M_A_X', 'ASR*', 'OLR*', 'ASR*-OLR*')
  
  %calculate the latitude of maximum energy flux
  flat=-89:.2:89;
  wid=25;
  mid=42;
  mn_lat=mid-wid;
  mx_lat=mid+wid;
  flat2=mn_lat:.1:mx_lat;
  for ny=1:NY;
      fAHT=interp1(mlat, squeeze(AHT(ny,:)), flat2)./cosd(flat2);
      AHT_max(ny,1) = mean(fAHT.*flat2)./mean(fAHT); %first moment of the heat transport weighted latitude over the range og latitudes
      AHT_max_mode(ny,1) = flat2(find2(cumsum(fAHT)/sum(fAHT), 0.5)); %mode of the AHT over the range 
      
      fAHT=interp1(mlat, squeeze(AHT(ny,:)), -flat2)./cosd(-flat2);
      AHT_max(ny,2) = mean(fAHT.*(-flat2))./mean(fAHT); %first moment of the heat transport weighted latitude over the range og latitudes
      AHT_max_mode(ny,2) = -flat2(find2(cumsum(fAHT)/sum(fAHT), 0.5)); %mode of the AHT over the range 
  end
  
 load('F:\trace\scripts\jet_location.mat','U_max', 'YEARS')
  figure(5);clf;
  hold on; plot(YEARS, -AHT_max(:,2), 'b', 'LineWidth', LW);
  hold on; plot(YEARS, -U_max(:,2), 'k', 'LineWidth', LW);
  hold on; plot(YD*[1 1], [ymn ymx], 'k--', 'LineWidth', 2)
  ymn=40; ymx=44;
  
  
  figure(6); hold on; 
  hold on; plot(YEARS, AHT_max(:,1), 'r', 'LineWidth', LW);
  hold on; plot(YEARS, U_max(:,1), 'k', 'LineWidth', LW);
  hold on; plot(YD*[1 1], [ymn ymx], 'k--', 'LineWidth', 2)



  
 
 
   