clear; 

% load obs parameter
load /data/observations/enso.reosc.parameter.seasonal.mat
a21o(1:12) = mean(a21);
a12o(1:12) = mean(a12);
t_obs = sst;
h_obs = tcline;




% numerical parameters
tdim = 2e2; % number of years
time = (1:12*tdim)/12;

% westerly wind events
wwe=zeros(1,12*tdim);

% initital values
T0 = 0.0;
h0 = 0.0;

%------------------- control run -------------------
% model parameter
a11c(1:12) = a11;
a21c(1:12) = a21;
a12c(1:12) = a12;
a22c(1:12) = a22;

scl     = 0.7;
noise_tc = scl*nonlinear1;
noise_hc = scl*nonlinear2;

nseed = 4; % pick a set of random numbers / set to fix value if always the same
%nseed = floor(1e4*rand); % each time pick a random set of random numbers 
seed=rng(nseed);         % set random numbers
[t_ctrl,h_ctrl]=reosc_model(T0,h0,a11c,a12c,a21c,a22c,noise_tc,noise_hc,wwe,tdim);


%------------------- sensitivity run -------------------
% model parameter
a11s(1:12) = mean(a11);
a21s(1:12) = mean(a21);
a12s(1:12) = mean(a12);
a22s(1:12) = mean(a22);

scl     = 0.7;
noise_ts(1:12) = scl*mean(nonlinear1);
noise_hs(1:12) = scl*mean(nonlinear2);

seed=rng(nseed);         % set random numbers
[t_sens,h_sens]=reosc_model(T0,h0,a11s,a12s,a21s,a22s,noise_ts,noise_hs,wwe,tdim);

% seasonal statistics
for m=1:12
    stdv_obs(m)  = std(t_obs(m:12:276));
    stdv_ctrl(m) = std(t_ctrl(m:12:12*tdim));
    stdv_sens(m) = std(t_sens(m:12:12*tdim));
end


% plot results

% initial circle
[tcc,hcc]=reosc_model(2,0,0.*(1:12),a12o,a21o,0.*(1:12),0.*(1:12),0.*(1:12),wwe,5);


% Control Experiment plots
figure(1); clf
subplot(3,1,1)
plot(10*a11c,'r'); hold on
plot(40*a12c,'b')
plot(a21c,'g')
plot(10*a22c,'c')
plot(20*noise_tc,'r--')
plot(0.2*noise_hc,'b--')
xlabel('calendar month'); ylabel('parameter values')
legend('a_{11} (growth rate T; x10)',...
       'a_{12} (coupling of T to h; x40)', ...
       'a_{21} (coupling of h to T; x10)', ...
       'a_{22} (growth rate h; x10)', ...
       'x1 (noise for T; x20)', ...
       'x2 (noise for h; x0.2)',-1 )
axis([1 12 -2 2]); grid on 
title('Model parameter seasonally')
text(7, 2.7,'Control Experiment','fontsize',15)

subplot(3,1,2)
plot(time,t_ctrl,'r'); hold on
plot(time,0.1*h_ctrl,'b')
xlabel('time [yrs]'); ylabel('T [^oC]     h [3m]')
legend('T','h')
axis([tdim-50 tdim -3 3]); grid on 
title('Last 50yrs of monthly mean T and h')

subplot(3,2,5)
plot(t_ctrl,h_ctrl); hold on
plot(tcc,hcc,'k','linewidth',2)
xlabel('T [^oC]'); ylabel('h [m]')
title('Phase diagram of monthly mean T and h')
axis([-4 4 -25 25]); grid on 

subplot(3,2,6)
plot(stdv_obs,'k'); hold on
plot(stdv_ctrl)
plot(stdv_sens,'r')
legend('observed','control','sense. exp.',4)
axis([1 12 0 1.5]); grid on
title('SST seasonal standard deviation')
xlabel('calendar month'); ylabel('T standard deviation [^oC]')


% Sensitivity Experiment plots
figure(2); clf
subplot(3,1,1)
plot(10*a11s,'r'); hold on
plot(40*a12s,'b')
plot(a21s,'g')
plot(10*a22s,'c')
plot(20*noise_ts,'r--')
plot(0.2*noise_hs,'b--')
xlabel('calendar month'); ylabel('parameter values')
legend('a_{11} (growth rate T; x10)',...
       'a_{12} (coupling of T to h; x40)', ...
       'a_{21} (coupling of h to T; x10)', ...
       'a_{22} (growth rate h; x10)', ...
       'x1 (noise for T; x20)', ...
       'x2 (noise for h; x0.2)',-1 )
axis([1 12 -2 2]); grid on 
title('Model parameter seasonally')
text(7, 2.7,'Sensitivity Experiment','fontsize',15)

subplot(3,1,2)
plot(time,t_sens,'r'); hold on
plot(time,0.1*h_sens,'b')
xlabel('time [yrs]'); ylabel('T [^oC]     h [3m]')
legend('T','h')
axis([tdim-50 tdim -3 3]); grid on 
title('Last 50yrs of monthly mean T and h')
text(9, 4,'Control Experiment','fontsize',15)

subplot(3,2,5)
plot(t_sens,h_sens); hold on
plot(tcc,hcc,'k','linewidth',2)
xlabel('T [^oC]'); ylabel('h [m]')
title('Phase diagram of monthly mean T and h')
axis([-4 4 -25 25]); grid on 

subplot(3,2,6)
plot(stdv_obs,'k'); hold on
plot(stdv_ctrl)
plot(stdv_sens,'r')
legend('observed','control','sense. exp.',4)
axis([1 12 0 1.5]); grid on
title('SST seasonal standard deviation')
xlabel('calendar month'); ylabel('T standard deviation [^oC]')


return
