%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function do_all=Stommel3box
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a Stommel-like 3box model, for APM 115 workshop
% Eli Tziperman
close all; clear all;

p=set_parameters;
ode_options=odeset('RelTol',1.e-7,'abstol',1.e-8);

% find steady state to initialize hysteresis run:
p.hysteresis=0;
tspan=[1 20000 21000]*p.year;
X0=[30 30 30]'; % (S1 S2 S3)
p.Temperature=[30 15 0];
[time,X] = ode45(@(t,X) rhs(t,X,p),tspan,X0,ode_options);
X0=X(3,:);
fprintf('initial q12=%g; q23=%g\n',q12(X0,p)/p.Sverdrup,q23(X0,p)/p.Sverdrup)

% do hysteresis run:
p.hysteresis=1;
tspan=[0:p.time_max/1000:p.time_max];
[time,X] = ode45(@(t,X) rhs(t,X,p),tspan,X0,ode_options);


% calculate all equilibria (stable and unstable) using fsolve:
p.hysteresis=1;
p.optimization=1;
fsolve_options=optimset('display','off','tolx',1.e-10,'tolfun',1.e-10);
ii=0;
for t=0:p.time_max/150:p.time_max;
  ii=ii+1;
  Fs=FW(t,p); FW_fsolve(ii)=Fs(1)*p.year/p.S0;
  for jj=1:p.Nrand
    X0=20*(rand(3,1)-0.5)+35.0;
    XX = fsolve(@(X)rhs(t,X,p),X0,fsolve_options);
    norm_b=norm(rhs(t,X0,p));
    norm_a=norm(rhs(t,XX,p));
    % fprintf('ii=%g; rhs before=%g; after=%g\n',ii,norm_b,norm_a)
    if norm_a>1.e-5
      Xfsolve(ii,jj,:)=NaN*XX(:);
    else
      Xfsolve(ii,jj,:)=XX(:)-mean(XX(:))+35.0;
    end
  end
end


% plots:
figure
subplot(3,4,1)
plot(time/p.year,X(:,1),'r')
hold on
plot(time/p.year,X(:,2),'g')
plot(time/p.year,X(:,3),'b')
% plot only for part of the simulation:
xlabel('time (yr)');
ylabel('S');

subplot(3,4,2)
i=0;
for t=tspan
  i=i+1;
  q12plot(i)=q12(X(i,:),p)/p.Sverdrup;
  q23plot(i)=q23(X(i,:),p)/p.Sverdrup;
end
plot(time/p.year,q12plot,'r')
hold on
plot(time/p.year,q23plot,'g')
% plot only for part of the simulation:
xlabel('time (yr)');
ylabel('q12,q23');

subplot(3,4,3)
i=0;
p.hysteresis=1;
for t=tspan
  i=i+1;
  Fs=FW(t,p);
  FWplot(i)=Fs(1)*p.year/p.S0;
end
plot(time/p.year,FWplot)
% plot only for part of the simulation:
xlabel('time (yr)');
ylabel('FW (m/yr)');

subplot(3,4,4)
N=length(q12plot); N2=floor(N/2);
plot(FWplot(1:N2),q12plot(1:N2),'r.')
hold on
plot(FWplot(N2+1:N),q12plot(N2+1:N),'b.')
% plot only for part of the simulation:
title('hysteresis');
xlabel('FW (m/yr)');
ylabel('q12 (Sv)');

subplot(3,4,5)
N=length(q23plot); N2=floor(N/2);
plot(FWplot(1:N2),q23plot(1:N2),'r.')
hold on
plot(FWplot(N2+1:N),q23plot(N2+1:N),'b.')
% plot only for part of the simulation:
title('hysteresis');
xlabel('FW (m/yr)');
ylabel('q23 (Sv)');

subplot(3,4,6)
N=length(X(:,1)); N2=floor(N/2);
plot(FWplot(1:N2),X(1:N2,1),'r.')
hold on
plot(FWplot(N2+1:N),X(N2+1:N,1),'b.')
% plot only for part of the simulation:
title('hysteresis');
xlabel('FW (m/yr)');
ylabel('S1');

subplot(3,4,7)
N=length(X(:,2)); N2=floor(N/2);
plot(FWplot(1:N2),X(1:N2,2),'r.')
hold on
plot(FWplot(N2+1:N),X(N2+1:N,2),'b.')
% plot only for part of the simulation:
title('hysteresis');
xlabel('FW (m/yr)');
ylabel('S2');

subplot(3,4,8)
N=length(X(:,3)); N2=floor(N/2);
plot(FWplot(1:N2),X(1:N2,3),'r.')
hold on
plot(FWplot(N2+1:N),X(N2+1:N,3),'b.')
% plot only for part of the simulation:
title('hysteresis');
xlabel('FW (m/yr)');
ylabel('S3');

subplot(3,4,9)
N=length(Xfsolve(:,1,1));
for ii=1:N;
  for jj=1:p.Nrand
    plot(FW_fsolve(ii),q12(Xfsolve(ii,jj,:),p)/p.Sverdrup,'r.');
    hold on
  end
end
% plot only for part of the simulation:
title('fsolve');
xlabel('FW (m/yr)');
ylabel('q12');

subplot(3,4,10)
N=length(Xfsolve(:,1,1));
for ii=1:N;
  for jj=1:p.Nrand
    plot(FW_fsolve(ii),q23(Xfsolve(ii,jj,:),p)/p.Sverdrup,'r.');
    hold on
  end
end
% plot only for part of the simulation:
title('fsolve');
xlabel('FW (m/yr)');
ylabel('q23');

figure
subplot(1,3,1)
N=length(Xfsolve(:,1,1));
for ii=1:N;
  for jj=1:p.Nrand
    plot(FW_fsolve(ii),Xfsolve(ii,jj,1),'r.');
    hold on
  end
end
% plot only for part of the simulation:
title('fsolve');
xlabel('FW (m/yr)');
ylabel('S1');

subplot(1,3,2)
N=length(Xfsolve(:,1,1));
for ii=1:N;
  for jj=1:p.Nrand
    plot(FW_fsolve(ii),Xfsolve(ii,jj,2),'r.');
    hold on
  end
end
% plot only for part of the simulation:
title('fsolve');
xlabel('FW (m/yr)');
ylabel('S2');

subplot(1,3,3)
N=length(Xfsolve(:,1,1));
for ii=1:N;
  for jj=1:p.Nrand
    plot(FW_fsolve(ii),Xfsolve(ii,jj,3),'r.');
    hold on
  end
end
% plot only for part of the simulation:
title('fsolve');
xlabel('FW (m/yr)');
ylabel('S3');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function flow=q12(X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DeltaT=p.Temperature(2)-p.Temperature(1);
DeltaS=X(2)-X(1);
flow=-p.k*(p.alpha*DeltaT-p.beta*DeltaS);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function flow=q23(X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DeltaT=p.Temperature(3)-p.Temperature(2);
DeltaS=X(3)-X(2);
flow=-p.k*(p.alpha*DeltaT-p.beta*DeltaS);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r=rhs(time,X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

deltaS12=X(2)-X(1);
deltaS23=X(3)-X(2);
Fs=FW(time,p);

r=[Fs(1)*p.area+abs(q12(X,p))*deltaS12, ...
   Fs(2)*p.area-abs(q12(X,p))*deltaS12+abs(q23(X,p))*deltaS23, ...
   Fs(3)*p.area-abs(q23(X,p))*deltaS23 ]'/p.volume;

if p.optimization==1
  r=r*1.e12;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Fs=FW(time,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% surface salt flux in northern box

if p.hysteresis
  if time<(p.time_max/2)
    ramp=time/(p.time_max/2);
  else
    ramp=1-(time-(p.time_max/2))/(p.time_max/2);
  end
  amp=p.FW_amp;
else
  ramp=1;
  amp=p.FW_amp0;
end

Fs(1)=   ramp*p.S0*amp*p.meter/p.year;
Fs(2)=   ramp*p.S0*amp*p.meter/p.year;
Fs(3)=-2*ramp*p.S0*amp*p.meter/p.year;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=set_parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p.meter=1;
p.sec=1;
p.km=1000*p.meter;
p.area=(4000*p.km)^2;
p.depth=4*p.km;
p.volume=p.area*p.depth;
p.year=365*24*3600;
p.time_max=50000*p.year;
p.alpha=0.0004; % K^-1
p.beta=0.001; % ppt^-1
p.Sverdrup=10^6*p.meter^3/p.sec;
p.k=10*p.Sverdrup/0.005; % Sv/(delta sigma)
p.S0=35.0;
p.FW_amp=1.3; % meter/year for hysteresis runs
p.FW_amp0=0.2; % meter/year for calculating steady state that
               % initializes the hysteresis.
p.Nrand=50;
p.optimization=0;
