%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function do_all=Stommel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stommel model, for APM 115
%% Eli Tziperman
close all; clear all;

disp([' Please see a version that also solves for the unstable ' ...
      'steady states in the workshop on climate.']);

DeltaT=15;
meter=1;
year=365*24*3600;
time_max=2000000*year;

% first, find steady state to initialize hysteresis run:
hysteresis=0;
tspan=[1 20000 21000]*year; X0=[38 32];
[T,X] = ode45(@(t,X) rhs_S(t,X,time_max,DeltaT,hysteresis),tspan,X0);
X0=X(3,:);

% multiply by these scales for plotting:
qscale=10*1.e10;fwscale=1/(35*0.0003*meter/year);

% next, do hysteresis run:
hysteresis=1;
tspan=[0:time_max/1000:time_max];
options=odeset('reltol',1.e-8);
[T,X] = ode45(@(t,X) rhs_S(t,X,time_max,DeltaT,hysteresis),tspan,X0,options);

subplot(3,1,1)
i=0;
hysteresis=1;
for t=tspan
  i=i+1;
  FWplot(i)=FW(t,time_max,hysteresis);
end
plot(T/year/1000,FWplot*fwscale,'LineWidth',2)
% plot only for part of the simulation:
h1=xlabel('time (kyr)');
h2=ylabel('FW ("meter/year")');
h3=title('fresh water forcing (FW) for hysteresis run');
set([h1,h2,h3,gca],'fontsize',13);

subplot(3,1,2)
i=0;
for t=tspan
  i=i+1;
  qplot(i)=q(DeltaT,X(i,:));
end
plot(T/year/1000,qplot*qscale,'LineWidth',2)
% plot only for part of the simulation:
h1=xlabel('time (kyr)');
h2=ylabel('THC ("Sv")');
h3=title('THC transport');
V=axis; axis([V(1) V(2) -5 20]);
set([h1,h2,h3,gca],'fontsize',13);

subplot(3,1,3)
N=length(qplot); N2=floor(N/2);
plot(FWplot(1:N2)*fwscale,qplot(1:N2)*qscale,'r.')
hold on
plot(FWplot(N2+1:N)*fwscale,qplot(N2+1:N)*qscale,'b.')
% plot only for part of the simulation:
h1=title('THC vs FW: red-1st half of run; blue: 2nd');
h2=xlabel('FW');
h3=ylabel('THC');
set([h1,h2,h3,gca],'fontsize',13);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function flow=q(DeltaT,X)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THC transport as function of temperature and salinity difference
% between the boxes
alpha=0.0004; % K^-1
beta=0.001; % ppt^-1
k=3e-8; % sec^-1
DeltaS=X(2)-X(1);
flow=k*(alpha*DeltaT-beta*DeltaS);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rhs=rhs_S(time,X,time_max,DeltaT,hysteresis)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% right hand side of the salinity equations
lambda=3e-11; % sec^-1
DeltaS=X(2)-X(1);
rhs=[FW(time,time_max,hysteresis)+abs(q(DeltaT,X))*DeltaS, ...
     -FW(time,time_max,hysteresis)-abs(q(DeltaT,X))*DeltaS]';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Fs=FW(time,time_max,hysteresis)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% surface Fresh Water (salt) flux in northern box

meter=1;
year=365*24*3600;
S0=35.0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify maximum and minimum of freshwater forcing range during the
% run.  These are given in relative units, so the values of both needs
% to be between 0 and 1.  A value of 0.5 corresponds to present-day
% climate:
FW_min=0;
FW_max=1; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if hysteresis
  if time<(time_max/2)
    FW_amplitude=FW_min+(FW_max-FW_min)*time/(time_max/2);
  else
    FW_amplitude=FW_max-(FW_max-FW_min)*(time-(time_max/2))/(time_max/2);
  end
else
  FW_amplitude=0.5;
end
Fs=-FW_amplitude*S0*0.0003*meter/year;
