clear; 

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

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

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

% initital values
T0 =  1.0;
h0 =  0.0;


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

% model parameter
a11m(1:12) = -0.08;
a12m(1:12) =  0.02;
a21m(1:12) = -1.11;
a22m(1:12) = -0.021;

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

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

% model parameter
a11m(1:12) = -0.08;
a12m(1:12) =  0.02;
a21m(1:12) = -1.11;
a22m(1:12) = -0.021;

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


% plot results

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

clf
subplot(2,2,1)
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 -3 3]); grid on 
title('Time series of monthly mean T and h')
text(9, 4,'Control Experiment','fontsize',15)

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

subplot(2,2,3)
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 -3 3]); grid on 
title('Time series of monthly mean T and h')
text(9, 4,'Sensitivity Experiment','fontsize',15)

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

return