clear;
alpha=[7 2 .15 .25 ]*1e-7; 		% in 1/s
phibar=[8 3 1 3.5];  % ozone mass mixing ratio in ppm

V=1;  % m/s
L=4e6; % m

m=1:4;
m1=[4 1:3];
m2=[3:4 1:2];
m3=[2:4 1];

R=V/L;
Rj=R+alpha;
Rjk=Rj(m).*Rj(m1);
Rjkl=Rj(m).*Rj(m1).*Rj(m2);
R1234=prod(Rj);
Const=(1-R^4/R1234);

phi=( alpha(m).*phibar(m)./Rj(m) ...
    + R*alpha(m1).*phibar(m1)./Rjk(m) ...
    + R.^2*alpha(m2).*phibar(m2)./Rjkl(m) ...
    + R.^3*alpha(m3).*phibar(m3)./R1234   )/Const;

dp=2500; % pa
g=9.8; %m/s2
Mol_wt=0.048; % kg/mol
Av=6.022e23; % molecules/mol
x=[-1 1]*45;
C=1e-6*dp/g/Mol_wt*Av*1e-4;  % converts to molecules/cm^2

% sum boxes in vertical columns to get column densities 
col_den0=C*[sum(phibar([1 4])) sum(phibar(2:3))];
col_den1=C*[sum(phi([1 4])) sum(phi(2:3))];

phi

disp(sprintf('Column Ozone in molecules/cm^2: Left (Boxes 1+4) and Right (Boxes 2+3) Columns\n',col_den0))
disp(sprintf('  for chemical equilibrium: %e %e  (solid line)\n',col_den0))
disp(sprintf('  with transport: %e %e  (dashed line)\n',col_den1))

plot(x,col_den0,x,col_den1,'--')
hold;plot(x,col_den0,x,col_den1,'--');hold
xlabel('latitude');  ylabel('Column Ozone - molecules cm^{-2}')



