%
%		This is dopipe1
%		It models the flow of North Atlantic Surface Water
%		down a long pipe. It may be regaded as a stream tube
%		which outcrops at the ocean surface and flows down
%		into the thermocline and around the gyre
%
%		Flow speed is .01 m/s (1 cm/s), mixing is 1000 m2/s
%		Length of tube is 12,000,000 m (12,000 km)
%		Delta-X is 50000 m (50 km)
%		Delta-T is 200,000 s (about 2.5 days)
%		We'll run it for 32 years (1950-1981)
%
%	first get surface tritium history (t in first col, trit
%						in second col)
%
    load natrsf.dat; t=natrsf(:,1); Cs=natrsf(:,2);
%
%		# of seconds in a year & other constants
%
    tyr=24*3600*365.25; lambda = 1/18/tyr;
    dx=50000; dt=200000; u=.01; K=1000;
%
%		# of steps in space and time
%
    L=12000000; nt=32*tyr/dt; nx=L/dx;
%
%		compute the weights (w's)
%
    Wm=u*dt/2/dx + K*dt/dx/dx;
    W0=1 - 2*K*dt/dx/dx;
    Wp=-u*dt/2/dx + K*dt/dx/dx;
    decay=dt*lambda;
%
%		initialize the data arrays
%		note that we need "old" arrays to hold the
%		previous time step.
%
    C=zeros(nx,1); Cold=C; H=zeros(nx,1); Hold=H;
    clg;				% clear graphics
    x=dx*[1:nx]/1.e6;			% distance in MegaMeters
%
    oldnow=0;			% flag to say we've plotted
%
%		the main calculation loop
%
for n=1:nt
    now=ceil(n*dt/tyr);			% the year number
    if now ~= oldnow 			% is it a new year?
	oldnow=now;			% update the year count
	subplot(2,1,1); p=plot(x,C,'r'); set(p,'LineWidth',2);hold on
	set(gca,'LineWidth',2,'FontSize',16);	
	ylabel('Tritium'); title(sprintf('Year = %d',1950+now-1));
	xlabel('Distance Down Pipe (Mm)');
	subplot(2,1,2); p=plot(x,H,'r'); set(p,'LineWidth',2);hold on
	set(gca,'LineWidth',2,'FontSize',16);	
	ylabel('Helium-3'); xlabel('Distance Down Pipe (Mm)');
	pause(1)		% to show the plot
	if now == 23			% do this for 1972
	    subplot(2,1,1); p=plot(x,C,'g'); set(p,'LineWidth',2);
	    subplot(2,1,2); p=plot(x,H,'g'); set(p,'LineWidth',2);
	    C72=C; H72=H;
	end		% end of inner if block
    end		% end of outer if block
%
% 		set incoming water to surface concentrations
%
    C(1)=Cs(now); H(1)=0;
%
%	do time step (staggered index trick!!!)
%
    C(2:nx-1) = (W0-decay)*Cold(2:nx-1) + Wm*Cold(1:nx-2) + Wp*Cold(3:nx);
    H(2:nx-1) = W0*Hold(2:nx-1) + Wm*Hold(1:nx-2) + Wp*Hold(3:nx) + decay*Cold(2:nx-1);
%
%	update old arrays
%
    Cold=C; Hold=H;

end		% end of main loop

% Finished, now replot the 1972 and the 1981 data

subplot(2,1,1); p=plot(x,C72,'g'); set(p,'LineWidth',2);
		p=plot(x,C,'b'); set(p,'LineWidth',2);
subplot(2,1,2); p=plot(x,H72,'g'); set(p,'LineWidth',2);
		p=plot(x,H,'b'); set(p,'LineWidth',2);