clear; 


% load obs parameter
load ./enso.reosc.parameter.seasonal.mat

tdim = 5; % number of years
time = (1:12*tdim)/12;

% model parameter
noise_t(1:12) = 0.;
noise_h(1:12) = 0.;

% westerly wind events
amp_wwe =   1e-9; % forcing on T and h

% initital values
T0 =   0.0;
h0 =   0.0;

%-----   control run ---------- 

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

% model parameter
a11m(1:12) = a11;
a21m(1:12) = a21;
a12m(1:12) = a12;
a22m(1:12) = a22;

[t_ctrl,h_ctrl]=reosc_model(T0,h0,a11m,a12m,a21m,a22m,noise_t,noise_h,wwe1,tdim);

%----- sensitivity run ---------- 

% westerly wind events
wwe2=zeros(1,12*tdim);
wwe2(1:12) = 1*amp_wwe;

% model parameter
a11m(1:12) = a11;
a21m(1:12) = a21;
a12m(1:12) = a12;
a22m(1:12) = a22;

[t_sens,h_sens]=reosc_model(T0,h0,a11m,a12m,a21m,a22m,noise_t,noise_h,wwe2,tdim);


% plot results

% reference circle
a11m(1:12)=0.; a12m(1:12) =  0.02; a21m(1:12) = -1.11; a22m(1:12) =0.; wwe=zeros(1,12*tdim);
[tcc,hcc]=reosc_model(2,0,a11m,a12m,a21m,a22m,noise_t,noise_h,wwe,tdim);

clf
subplot(2,3,1)
plot(12*time,wwe1/amp_wwe,'bo'); hold on
xlabel('time [mon]'); ylabel('wind bursts [norm.]')
axis([1 12 -12 12]); grid on 
title('Wind burst in first 12 mon.')

subplot(2,3,2)
plot(time,t_ctrl,'r'); hold on
plot(time,0.3*h_ctrl,'b')
xlabel('time [yrs]'); ylabel('T [^oC]     h [3m]')
legend('T','h')
axis([0 tdim -4 4]); grid on 
title('Time series of monthly mean T and h')
text(0.5, 5,'Control Experiment','fontsize',15)

subplot(2,3,3)
plot(t_ctrl(1),h_ctrl(1),'o'); hold on
plot(t_ctrl(12:12:12*tdim),h_ctrl(12:12:12*tdim),'ro')
plot(t_ctrl(1),h_ctrl(1),'o')
plot(t_ctrl,h_ctrl)
plot(tcc,hcc,'k')
xlabel('T [^oC]'); ylabel('h [m]')
title('Phase diagram of monthly mean T and h')
legend('initital condition [1 month]','Dec',4)
axis([-4 4 -25 25]); grid on 

subplot(2,3,4)
plot(12*time,wwe2/amp_wwe,'bo'); hold on
xlabel('time [mon]'); ylabel('wind bursts [norm.]')
axis([1 12 -12 12]); grid on 
title('Wind burst in first 12 mon.')

subplot(2,3,5)
plot(time,t_sens,'r'); hold on
plot(time,0.3*h_sens,'b')
xlabel('time [yrs]'); ylabel('T [^oC]     h [3m]')
legend('T','h')
axis([0 tdim -4 4]); grid on 
title('Time series of monthly mean T and h')
text(0.5, 5,'Sensitivity Experiment','fontsize',15)

subplot(2,3,6)
plot(t_sens(1),h_sens(1),'o'); hold on
plot(t_sens(12:12:12*tdim),h_sens(12:12:12*tdim),'ro')
plot(t_sens(1),h_sens(1),'o')
plot(t_sens,h_sens)
plot(tcc,hcc,'k')
xlabel('T [^oC]'); ylabel('h [m]')
title('Phase diagram of monthly mean T and h')
legend('initital condition [1 month]')
axis([-4 4 -25 25]); grid on 

return