% Stochastic Resonance of the 2-box Stommel model driven by periodic
% plus stochastic FW forcing.

% Kate Dennis & Amanda Hughes (AM115 Project 2, spring 2007)
% Time given in years

%Achieve reversal with ramping up noise or amplitude of Milankovitch

function solvestommelsec % (span,initial)

span=[0 10000];
initial=[35 35];
[t,S]=ode45(@salinity,span,initial);

figure()
%Plot salinity in box 1 & 2 versus time
subplot(2,1,1)
plot(t,S(:,1),'r','LineWidth',1)
hold on
plot(t,S(:,2),'b','LineWidth',1)
xlabel('Time x 10^3 years')
ylabel('Salinity')
legend('Low latitude salinity', 'High latitude salinity')

%Plot flux versus time
subplot(2,1,2)

T(1) = 303; %degrees K
T(2) = 278; %degrees K
k=50;%*Sv/0.005; %m3 per second [Sverdrup (x10^6 m3 s-1)/delta rho]
alpha = 0.0004; %per K
beta = 0.001; %per ppt

DeltaT = T(2)-T(1);
DeltaS = S(:,2)-S(:,1);
flow = k*(-alpha*DeltaT+beta*DeltaS); %m3 per second

plot(t,flow,'LineWidth',1)
xlabel('Time x 10^3 years')
ylabel('Flux (Sverdrup)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dS = salinity(t,S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%yr = 365*24*3600; %Number of seconds in one year
%A=1.6^13; %Surface area of one box (m2)
%volume=A*4000; %Volume of one box (m3)
Fs = 1;%*A*yr/volume; %m3 per year
Flux12 = Q12(S); %Density driven flux between box 1 and box 2 [Sverdrup(x10^6m3 per second)]


% Change these two parameters for amplitude of periodic and
% stochastic forcing components:
% ----------------------------------------------------------
% Mamp=Amplitude of periodic forcing
% Samp=Amplitude of stochasicity

% A few interesting cases:
% no transitions:
% Mamp = 0.5; Samp = 3;
% SRa) stocahstic resonance: some transitions, some skipped periods:
 Mamp = 0.5; Samp = 4;
% SRb) more transitions, still weakly periodic, stochastic transitions:
% Mamp = 0.5; Samp = 6;
% SRc) stochastic transitions:
% Mamp = 0.5; Samp = 12;

% ----------------------------------------------------------

rand('state',sum(1000*round(t)));
F=(Fs-Samp/2+Samp*rand()+Mamp*cos(2*pi/1000*t));%salinity*m3 per year  
dS=zeros(2,1);

dS(1) = F+abs(Flux12)*(S(2)-S(1));   
dS(2) = -F+abs(Flux12)*(S(1)-S(2)); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function flow = Q12(S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%yr = 365*24*3600; %Number of seconds in one year
%Sv = 10^6; %Convert k from Sverdrups to SI units (x10^6 m3 s-1)
%A = 1.6^13; %Surface area of one box (m2)
%volume=A*4000; %Volume of one box (m3)
T(1) = 303; %degrees K
T(2) = 278; %degrees K
k=100;%*Sv*yr/volume;%[Sverdrup (x10^6 m3 s-1)/delta rho]
alpha = 0.0004; %per K
beta = 0.001; %per ppt

DeltaT = T(2)-T(1);
DeltaS = S(2)-S(1);
flow = k*(-alpha*DeltaT+beta*DeltaS); %m3 per year


