%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
 
 nni=find(mlat == 0);
 AHT_EQ=AHT(:,nni);
 OHT_EQ=OHT(:,nni);
 RAD_EQ=OHT_EQ+AHT_EQ;
 RAD_EQ2=SW_HT(:,nni)+LW_HT(:,nni);
 cut_per=200/10;
 
  [low_AHT]=low_filter(AHT_EQ,1,cut_per); %200 year low pass filter
  [low_OHT]=low_filter(OHT_EQ,1,cut_per);
  [low_RAD]=low_filter(RAD_EQ,1,cut_per);
  [low_RAD2]=low_filter(RAD_EQ2,1,cut_per);
 
 LWd=1.3; LW=2.5; LWM=2.3;
 FSz=22; FN='Arial';
 PW=10^-15;
 %colors
lgreen=[.7 1 .7];
dgreen=[34 139 34]/255; 
lorange = [255 165 0]/255;
dorange=[255 100 0]/255;
purple=[160 32 240]/255;
lblue = [135 206 250]/255
 
 figure(1); clf;
 hold on; plot(YEARS, PW*AHT_EQ, 'b', 'LineWidth', LWd);
 hold on; plot(YEARS, PW*RAD_EQ, 'Color',lorange, 'LineWidth', LWd);
 hold on; plot(YEARS, -PW*OHT_EQ, 'Color',lgreen, 'LineWidth', LWd);
 
 hold on; plot(YEARS, PW*low_AHT, 'b', 'LineWidth', LW);
 hold on; plot(YEARS, PW*low_RAD, 'Color',dorange, 'LineWidth', LW);
 hold on; plot(YEARS, -PW*low_OHT, 'Color',dgreen, 'LineWidth', LW);
 
 hold on; plot(YEARS, PW*mean(low_AHT)*ones(size(YEARS)), 'b--', 'LineWidth', LWM);
 hold on; plot(YEARS, PW*mean(low_RAD)*ones(size(YEARS)), '--','Color',dorange, 'LineWidth', LWM);
 hold on; plot(YEARS, -PW*mean(low_OHT)*ones(size(YEARS)), '--','Color',dgreen, 'LineWidth', LWM);
 xticl=-20:5:0;
 xtic=xticl*1000;
 xlabel('Years X1000', 'FontSize', FSz, 'FontName', FN);
 ylabel('AHT_E_Q and hemispheric contrsat in eflux (PW)', 'FontSize',FSz, 'FontName', FN);
 set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2.5, 'XTick', xtic, 'XTickLabel', xticl, 'YTickLabel', []);
 box on;
 
 do=1:NY;
 CSHF=squeeze(mean(SHF(do,:,:)));
 CTOA=squeeze(mean(TOA(do,:,:)));
 COLR=squeeze(mean(OLR(do,:,:)));
 CASR=squeeze(mean(FSNT(do,:,:)));
 
 for ny=1:length(do);
     aSHF(ny, 1:NN, 1:N)=squeeze(SHF(do(ny), :,:)) - CSHF;
     aTOA(ny, 1:NN, 1:N)=squeeze(TOA(do(ny), :,:)) - CTOA;
     aASR(ny, 1:NN, 1:N)=squeeze(FSNT(do(ny), :,:)) - CASR;
     aOLR(ny, 1:NN, 1:N)=squeeze(OLR(do(ny), :,:)) - COLR;
 end

 ts=low_OHT;
 ts=(ts-mean(ts))/std(ts);
 for nn=1:NN;
     for n=1:N;
         SHF_MAP(nn,n) = mean(squeeze(aSHF(:,nn,n)).*ts);
     end
 end

CMAX=60;
load('F:\trace\LH_dec.mat', 'lat', 'lon') 
N_std=2;
sphcolor2(-N_std*double(SHF_MAP),lat,lon,4); %shading interp;
hold on; plotm(0*lon,lon,'m-', 'LineWidth', 2.3);
load coast;
hold on; plotm(lat,long,'k-', 'LineWidth', 2.3);
box on;  framem('on');tightmap2;
caxis([-CMAX CMAX]);
H=colorbar('peer', gca, 'LineWidth', 2.5);
set(get(H,'YLabel'),'string', char('Surface Heat Flux (W m^-^2)'), 'FontSize', FSz, 'FontName', FN);
load('F:\trace\LH_dec.mat', 'lat', 'lon') 

 %ocean basin masks
MASK=zeros(NN,N);
sphcolor2jet(MASK,lat,lon,3);
minlon=270*ones(NN,1);
do1=find(lat<50 & lat>=22);
minlon(do1) = 270 + 0.5*(lat(do1)-50); %coming into panama
do2=find(lat<=22 & lat>=0);
minlon(do2)=258-1.8*(lat(do2)-22);

do3=find(lat<=0 &lat >=-35);
minlon(do3)=296;

do4=find(lat<=-35 &lat >=-55);
minlon(do4)=296 +0.70*(lat(do4)+35);

do5=find(lat<=-55);
minlon(do5)= 283;

hold on; plotm(lat,minlon, 'k', 'LineWidth', 2);


maxlon=0*ones(NN,1);
do1=find(lat<=40 & lat>=-40);
maxlon(do1)=0 - 0.35*(lat(do1)-40);

do2=find(lat<=-40);
maxlon(do2)=28;

hold on; plotm(lat,maxlon, 'k', 'LineWidth', 2);

%now add in the MASK in the Atlantic by stepping between the min and max
%longitude
for nn=1:NN;
    do=find(lon>=minlon(nn) | lon<=maxlon(nn));
    AMASK(nn,do)=1;
end

sphcolor2jet(AMASK,lat,lon,4);

do1=find(lat>=20);
pmin_lon(do1)= 110*ones(size(do1));
figure(4); hold on;
plotm(lat(do1), pmin_lon(do1)', 'm-', 'LineWidth', 2)

do2=find(lat<=20);
pmin_lon(do2)=110-(50/90)*(lat(do2)-20);
plotm(lat, pmin_lon', 'm-', 'LineWidth', 2)
PMASK=zeros(NN,N);
IMASK=zeros(NN,N);
 for nn=1:NN;
     do=find(lon>=pmin_lon(nn) & lon<=minlon(nn));
     PMASK(nn,do)=1;
     
     do2=find(lon>=maxlon(nn) & lon<=pmin_lon(nn));
     IMASK(nn,do2)=1;
 end
 

 % sphcolor2jet(PMASK,lat,lon,5);

  %sphcolor2jet(IMASK,lat,lon,6);


%now zonally integrate and back out the ocean heat transport anomalies
G_SHF = mean(SHF_MAP');
A_SHF= mean((SHF_MAP.*AMASK)');
P_SHF= mean((SHF_MAP.*PMASK)');
I_SHF= mean((SHF_MAP.*IMASK)');
norm=trapz(sind(lat),ones(size(lat)));
GSHF=trapz(sind(lat),G_SHF)/norm;
ASHF=trapz(sind(lat),A_SHF)/norm;
ISHF=trapz(sind(lat),I_SHF)/norm;
PSHF=trapz(sind(lat),P_SHF)/norm;

GOHT(1)=0; AOHT(1)=0; IOHT(1)=0; POHT(1)=0;
for nn=2:length(mlat);
GOHT(nn) = GOHT(nn-1) + dA(nn-1)*(-G_SHF(nn-1)+GSHF);
AOHT(nn) = AOHT(nn-1) + dA(nn-1)*(-A_SHF(nn-1)+ASHF);
IOHT(nn) = IOHT(nn-1) + dA(nn-1)*(-I_SHF(nn-1)+ISHF);
POHT(nn) = POHT(nn-1) + dA(nn-1)*(-P_SHF(nn-1)+PSHF);
end

figure(7);clf; subplot(1,2,1);
SF=-N_std*10^-15;
hold on; plot(SF*GOHT,mlat, 'k', 'LineWidth', LW);
hold on; plot(SF*AOHT,mlat, 'r', 'LineWidth', LW);
hold on; plot(SF*POHT,mlat, 'b', 'LineWidth', LW);
hold on; plot(SF*IOHT, mlat, 'g', 'LineWidth', LW);
hold on; plot(zeros(size(mlat)),mlat, 'm-', 'LineWidth', LW);
hold on; plot([-0.25 0.1], [0 0], 'm', 'LineWidth', LW);
legend('Global', 'Atlantic', 'Pacific', 'Indian')
xlabel('Ocean Heat Transport Anomaly (PW)', 'FontSize', FSz, 'FontName', FN);
ylabel('Latitude', 'FontSize', FSz, 'FontName', FN);
set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2.2);
axis([-0.25 0.1 -90 90]); box on;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RAD maps
 ts=low_RAD;
 ts=(ts-mean(ts))/std(ts);
 for nn=1:NN;
     for n=1:N;
         RAD_MAP(nn,n) = mean(squeeze(aTOA(:,nn,n)).*ts);
         ASR_MAP(nn,n) = mean(squeeze(aASR(:,nn,n)).*ts);
         OLR_MAP(nn,n) = mean(squeeze(aOLR(:,nn,n)).*ts);
     end
 end

 [LON LAT]=meshgrid(lon,lat);
 LWC=2.4;
sphcolor2(N_std*double(RAD_MAP),lat,lon,7); %shading interp;
hold on; contourm(LAT,LON,N_std*double(ASR_MAP), [-100:20:-20], 'b', 'LineWidth', LWC);
hold on; contourm(LAT,LON,N_std*double(ASR_MAP), [20:20:100], 'r', 'LineWidth', LWC);
hold on; plotm(0*lon,lon,'m-', 'LineWidth', 2.3);
load coast;
hold on; plotm(lat,long,'k-', 'LineWidth', 2.3);
box on;  framem('on');tightmap2;
caxis([-CMAX CMAX]);
H=colorbar('peer', gca, 'LineWidth', 2.5);
set(get(H,'YLabel'),'string', char('TOA rad anomaly (W m^-^2)'), 'FontSize', FSz, 'FontName', FN);
load('F:\trace\LH_dec.mat', 'lat', 'lon') 


% sphcolor2(N_std*double(OLR_MAP),lat,lon,8); %shading interp;
% hold on; plotm(0*lon,lon,'m-', 'LineWidth', 2.3);
% load coast;
% hold on; plotm(lat,long,'k-', 'LineWidth', 2.3);
% box on;  framem('on');tightmap2;
% caxis([-CMAX CMAX]);
% H=colorbar('peer', gca, 'LineWidth', 2.5);
% set(get(H,'YLabel'),'string', char('OLR anomaly(W m^-^2)'), 'FontSize', FSz, 'FontName', FN);
% load('F:\trace\LH_dec.mat', 'lat', 'lon') 
% 
% sphcolor2(N_std*double(ASR_MAP),lat,lon,9); %shading interp;
% hold on; plotm(0*lon,lon,'m-', 'LineWidth', 2.3);
% load coast;
% hold on; plotm(lat,long,'k-', 'LineWidth', 2.3);
% box on;  framem('on');tightmap2;
% caxis([-CMAX CMAX]);
% H=colorbar('peer', gca, 'LineWidth', 2.5);
% set(get(H,'YLabel'),'string', char('ASR anomaly(W m^-^2)'), 'FontSize', FSz, 'FontName', FN);
% load('F:\trace\LH_dec.mat', 'lat', 'lon') 

%zonal mean radiation and shf
MAX=8.5;
figure(30);clf;subplot(1,2,1);
LW=2.7;
LW2=1.8;
hold on; plot(-N_std*double(mean(SHF_MAP')),lat, 'Color', dgreen,'LineWidth',LW ); %shading interp;
hold on; plot(N_std*double(mean(RAD_MAP')),lat, 'Color', dorange,'LineWidth',LW ); %shading interp;
%hold on; plot(N_std*double(mean(ASR_MAP')),lat, '--','Color','r' ,'LineWidth',LW ); %shading interp;
%hold on; plot(-N_std*double(mean(OLR_MAP')),lat, '--','Color','b' ,'LineWidth',LW ); %shading interp;

hold on; plot([-MAX MAX], [0 0], 'm', 'LineWidth', 2.3)
hold on; plot([0 0], [-90 90], 'm', 'LineWidth', 2.3)

axis([-MAX MAX -90 90]);

set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth', 2.3, 'YTick', [-80:20:80])
box on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


break

 %first look for the spikes over the late holocene
 do=find(YEARS >= -5000 & YEARS <= -500);
 CSHF=squeeze(mean(SHF(do,:,:)));
 CTOA=squeeze(mean(TOA(do,:,:)));
 COLR=squeeze(mean(OLR(do,:,:)));
 CASR=squeeze(mean(FSNT(do,:,:)));
 
 for ny=1:length(do);
     aSHF(ny, 1:NN, 1:N)=squeeze(SHF(do(ny), :,:)) - CSHF;
     aTOA(ny, 1:NN, 1:N)=squeeze(TOA(do(ny), :,:)) - CTOA;
     aASR(ny, 1:NN, 1:N)=squeeze(FSNT(do(ny), :,:)) - CASR;
     aOLR(ny, 1:NN, 1:N)=squeeze(OLR(do(ny), :,:)) - COLR;
 end
 
 ts=-OHT_EQ(do); ts=(ts-mean(ts));%/std(ts);
 figure(2);clf;
 hold on; plot(YEARS(do), ts, 'k-', 'LineWidth', LW)
 
 nti=find (ts == min(ts));
 MAX=20;
sphcolor2(double(aSHF(nti,:,:)),lat,lon,4); title('SHF anomaly', 'FontSize', FSz, 'FontName', FN); caxis([-MAX MAX]);colorbar; 
sphcolor2(double(aTOA(nti,:,:)),lat,lon,5); title('TOA RADIATION', 'FontSize', FSz, 'FontName', FN); caxis([-MAX MAX]);colorbar; 
 
 
 