% Sylvia G. Dee
% SEA Analysis: PHYDA, monthly NINO34 reconstructions

% Purpose: Does the latest Palmyra data support the hypothesis that El Nino
% events happen more often than usual after strong eruptions?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%clc
clear all

% Load Datasets

%load JEG_graphics
load('Toohey.mat')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('phyda_ninoX_ensemble_05-Jan-2018_14:24:35.mat')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load all ensemble data
time=mon_r;
NINO34_ens=xa_nino34_ens;
%NINO12_ens=xa_nino12_ens;
%NINO3_ens=xa_nino3_ens;
%NINO4_ens=xa_nino4_ens;
%NINOgrad_ens=xa_ninodelT_ens;


%NJS: Calculate annual cycle from each ensemble member then remove it
n34_mon12=reshape(NINO34_ens,[12,2000,998]);
n34_nsc=bsxfun(@minus,n34_mon12,mean(n34_mon12,2));% remove climatology


%Calculate ensemble mean for each variable & save as annual averages

NINO34_yr=mean(reshape(NINO34_ens,[12,2000, 998]),3);
%NINO12_yr=mean(reshape(NINO12_ens,[12,2000, 998]),3);
%NINO3_yr=mean(reshape(NINO3_ens,[12,2000, 998]),3);
%NINO4_yr=mean(reshape(NINO4_ens,[12,2000, 998]),3);
%NINOgrad_yr=mean(reshape(NINOgrad_ens,[12,2000, 998]),3);

% DOESNT WORK
% NINO34_error=quantile(reshape(NINO34_ens,[12,2000,998]),[.05, .95],3); %
% Cd5err=NINO34_error(:,:,1);
% Cd95err=NINO34_error(:,:,2);

NINO34_error=reshape(NINO34_ens,[12,2000,998]);

% set vars for calculation

years=(1:1:2000);

%% Volcanic Forcing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

VtT=Toohey(:,1);
VdT=Toohey(:,2);

S_time=VtT;
Sb=VdT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find intervals where coral data doesn't exist:

Vt = intersect(years, S_time);           % Use intersection of Coral Real-Valued Datapoints in Time and Volc_Time
[V_Index,LOC] = ismember(Vt,S_time);   % LOC returns place where T-Volc is equal to Vt. Pull that data out.
Vd=Sb(LOC);                      % Corrected Volcanic Dataset--Matches places where coral has real-valued data.

% 1. Select Key dates: 

Erup_Forcing = 0.22;
key_dates=find(Vd>Erup_Forcing);

% Define terup, the array of years in which eruptions occurred.

terup = Vt(key_dates); 

%% Remove Seasonal Cycle from all data to highlight yearly response to volc

% % ensemble mean
% seasons=mean(NINO34_yr,2);
% sea_array=repmat(seasons,[1,2000]);
% 
% nino34_m=NINO34_yr-sea_array;
% 
% % full ensemble -- remove seasonal cycle from each member individually
% nino34_ens_full=zeros(size(NINO34_error));
% 
% for i=1:998
%     seasons_err=mean(NINO34_error(:,:,i),2);
%     sea_arr=repmat(seasons_err,[1,2000]);
%     nino34_ens_full(:,:,i)=NINO34_error(:,:,i)-sea_arr;
% end

%% Plot Monthly values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(2,1,1)

%for i=1:998
%   plot(1:1:12,mean(nino34_ens_full(:,key_dates,i),2),'Color',ligr,'LineWidth',1.0),hold on;
%end

% Plot each of the ensemble members
for i=1:998;plot(n34_nsc(:,key_dates,i),'color',[0.4 0.4 0.4]);hold on; end
% Plot the ensemble mean for each eruption
n34_ensmn=mean(n34_nsc(:,key_dates,:),3);
plot(1:1:12,n34_ensmn,'r','LineWidth',3.0)
% Plot the ensemble mean, volcanic mean
plot(1:1:12,mean(n34_ensmn,2),'k','LineWidth',3.0)


%plot(1:1:12,mean(nino34_m(:,key_dates),2),'r','LineWidth',3.0),hold on;
%plot(1:1:12,mean(nino34_m,2), 'k','LineWidth',3.0)

title('A. YR 0 Volcanic Response, Monthly NINO34 SST, PHYDA')
ylabel('NINO3.4 SST ANOM')
xlabel('YEAR 0')

xticks([1:2:12])
xticklabels({'APR','JUN','AUG','OCT','DEC','FEB'})

x=[0,12];
grid on
y=[0,0];
set(gca,'FontName','Helvetica Neue','fontsize',18)
plot(x,y,'k--')
%legend('Eruption Years (AOD > 0.22)','Mean, All Years','location','southeast')
ylim([-2,4]);
xlim([1 12])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,1,2)
%for i=1:length(key_dates)

% for i=1:998
%     plot(1:1:12,mean(nino34_ens_full(:,key_dates+1,i),2),'Color',ligr,'LineWidth',1.0),hold on;
% end
% 
% plot(1:1:12,mean(nino34_m(:,key_dates+1),2),'r','LineWidth',3.0),hold on;
% plot(1:1:12,mean(nino34_m,2), 'k','LineWidth',3.0)


% Plot each of the ensemble members
for i=1:998;plot(n34_nsc(:,key_dates+1,i),'color',[0.4 0.4 0.4]);hold on; end
% Plot the ensemble mean for each eruption
n34_ensmn=mean(n34_nsc(:,key_dates+1,:),3);
plot(1:1:12,n34_ensmn,'r','LineWidth',3.0)
% Plot the ensemble mean, volcanic mean
plot(1:1:12,mean(n34_ensmn,2),'k','LineWidth',3.0)


xticks([1:2:12])
xticklabels({'APR','JUN','AUG','OCT','DEC','FEB'})


ylim([-2,4]);
xlim([1 12])

x=[0,12];
y=[0,0];
grid on
plot(x,y,'k--')
title('B. YR+1 Volcanic Response, Monthly NINO34 SST, PHYDA')
ylabel('NINO3.4 SST ANOM')
%legend('Eruption Years (AOD > 0.22)','Mean, All Years','location','southeast')

set(gca,'FontName','Helvetica Neue','fontsize',18)
xlabel('YEAR +1')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

