% Stommel's 1961 model of convection in coupled boxes.
%
% Compute the temperature and salinity difference, here represented
% by y and x, between two well-stirred boxes that are both
% in contact with a reservoir that has (y, x) = (1, 1). Both
% boxes conduct y and x at rates 1 and delta (delta < 1), and the
% density difference between the boxes is given by d = -y + R*x, where
% we are assuming R > 1, typically. There can be advection (or
% flushing) between the boxes at a rate d*lambda, where the
% flushing does not depend upon the sign of the density anomaly.
%
% This code also has 1) flushing with inertia (not especially
% interesting), 2) a flickering or random temperature
% perturbation (slightly interesting), and 3) an oscillating
% reservoir temperature (hoping for but not yet finding chaos;
% may be present in parameter regimes not checked).
%
% This code has not been fully tested, but seems to reproduce
% S61's two cases fairly well when all three of the 'new' things
% are turned off, of course.
%
% Written by Jim Price, April 28, 1999. Public domain.
%
set(0,'DefaultLineLineWidth',1.2)
set(0,'DefaultTextFontSize',14)
set(0,'DefaultAxesLineWidth',1.6)
set(0,'DefaultAxesFontSize',14)
clear; format compact;
nn = 0;
% set the model parameters (follows S61, aside from 'new')
R = 2.0; % abs of the ratio of the expansion coefficients, x/y
delta = 1/6; % conduction rate of salinity wrt temperature
lambda = 0.2; % inverse non-d flushing rate; = inf for no flushing
q = 0.; % initial flushing rate (0 to 1) 'new'
qdelta = 100.; % time constant (inertia) for flushing; 'new'
% set = 1/dtau for equilibrium flushing as in S61
% set = 0.2 for slowly responding flushing
yres = 1.; % steady reservoir y, = 1 for S61 case 'new'
resosc = 0.; % amplitude of reservoir y oscillation 'new'
dtau = 0.01; % the time step of non-d time; 0.01 seems OK
nstep = 1500; % number of time steps; 1500 is usually
% enough to insure convergence to a steady state
%
% This model version is set up to do integration over a range
% of initial T,S or y,x from 0 to 1. The increment of T and S
% are set by delT and delS = 1/ni, where ni is the number of
% integrations (n1 = 1 to 20 is reasonable).
yres0 = yres;
ni = 6;
delT = 1/ni; % make ni = 1 or 2 to reduce the number of
% integrations to be done
delS = delT;
for n1=0:delT:1 % n1 and n2 are used to set the initial T,S
for n2=0:delS:1
if n1==0 | n1==1 | n2==0 | n2==1 % skip all non-boundary initial points
x(1) = n1; % set the initial temperature
y(1) = n2; % set the initial salinity
for m = 2:nstep
tau(m) = m*dtau; % the non-d time
% evaluate the reservoir temperature (y); note that
% this temperature is steady if resosc = 0. (the S61 case)
yres = yres0 + resosc*sin(tau(m)*pi);
% the first part of a second order R-K method; time step forward
% by half a time step to make a first guess at the new times
dr = abs(R*x(m-1) - y(m-1)); % the density anomaly
qequil = dr/lambda; % the equilibrium flushing
yh = y(m-1) + dtau*(yres - y(m-1))/2 - ...
dtau*y(m-1)*q/2; % time step the temperature
xh = x(m-1) + dtau*delta*(1 - x(m-1))/2 - ...
dtau*x(m-1)*q/2; % time step the salinity
qh = q + dtau*qdelta*(qequil - q)/2; % time step the flushing
% the second part; use the half time step values to make a full step
dr = abs(R*xh - yh);
qequil = dr/lambda;
y(m) = y(m-1) + dtau*(yres - yh) - ...
dtau*qh*yh;
x(m) = x(m-1) + dtau*delta*(1 - xh) - ...
dtau*qh*xh;
q = q + dtau*qdelta*(qequil - qh);
% now add on a flickering temperature if you want to (or comment out)
% tflickamp = 0.01; % set the amplitude here
% tflick = tflickamp*unifrnd(-1., 1.);
% y(m) = y(m) + tflick;
end % end of time step loop
d = R*x - y; % evaluate the density
if nn == 0;
nn = 1;
% make a time series plot of the first case only
figure(1)
clf reset
subplot(2,1,1)
plot(tau, x,'-', tau, y,'--')
legend('salinity', 'temperature')
ylabel('T, S diff, non-d')
title('Experiment 1,1')
subplot(2,1,2)
plot(tau, d)
xlabel('time, non-d')
ylabel('density diff')
% contour the density (or flushing rate), and add the (T,S)
% trajectories on top of the contours
ym = 0:0.1:1;
xm = 0:0.1:1;
for k1 = 1:11
for k2 = 1:11
dm(k1,k2) = (1/lambda)*(R*xm(k2) - ym(k1));
end
end
figure(2)
clf reset
dc = -10:2:20;
c = contour(xm, ym, dm, dc,'k');
clabel(c);
xlabel('salinity diff, non-d')
ylabel('temp diff, non-d')
hold on
end % the if on nn = 0
[m1 m2] = size(x);
% plot the individual trajectories.
% color code according to which equilibrium point
% the trajectory ends up on. this will likely have be
% reset if the model parameters (R, delta, lambda) are changed.
if d(m2) >= 0
plot(x, y, 'r--')
plot(x(m2), y(m2), '*r')
else
plot(x, y, 'g')
plot(x(m2), y(m2), '*g')
end
clear x y d
end % or to gate out non-boundary values of n1, n2
end % loop on n1
end % loop on n2
% make some plots to show where roots (equil. points) are
for k=1:60
f(k) = (k-30)*0.1;
lhs(k) = lambda*f(k);
rhs(k) = (R/(1 + abs(f(k))/delta)) - 1/(1 + abs(f(k)));
end
figure(3)
clf reset
plot(f, rhs, f, lhs)
xlabel('f, flow rate')
ylabel('lhs(f), rhs(f)')
title('roots of S61 model')
grid
% end of the script