%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function do_all=Stommel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stommel model under global warming scenario
% Eli Tziperman
close all; clear all;
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);
%ScreenSize is a four-element vector: [left, bottom, width, height]:
scrsz = get(0,'ScreenSize');
% position= [left, bottom, width, height]
figure('Position',[100 50 scrsz(3)/2 scrsz(4)-150])
subplot(2,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(2,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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.5;
FW_max=0.6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;