Program diffuse1;
        { define variables...}
const  Ncel = 10; {90}
var   i, j, Nmix                     : integer;
      Totx, Tott, Delx, Delt, De, mixf  : real;
      Ct1, Ct2                       : array[0..Ncel] of real;

begin
        { set variables and initial concentrations...}
  De := 1e-5 {cm2/s}; Totx := 1000.0 {cm}; Tott := 100 * 365 * 24 * 60 * 60.0 {s}; Delx := Totx / Ncel;
  for i := 1 to Ncel do Ct1[i] := 0.0; Ct1[0] := 0.56;
	{ calculate the mixing factor ...}
  mixf := De * Tott / (Delx*Delx);
  Nmix := 1 + trunc(mixf / 0.22);                    {trunc(..) gives integer part }
{  if Nmix < 3000 then Nmix := 3000;}
  mixf := mixf / Nmix;
        { now diffuse...}
  for i:= 1 to Nmix do
    begin                    { mix into Ct2-cells ...}
      Ct2[1] := 2 * mixf * Ct1[0] + mixf * Ct1[2] + (1 - 3 * mixf) * Ct1[1];         { .. first cells }
      for j := 2 to Ncel - 1 do
        Ct2[j] := mixf * (Ct1[j-1] + Ct1[j+1]) + (1 - 2 * mixf) * Ct1[j];         { .. inner cells }
      Ct2[Ncel] := mixf * Ct1[Ncel-1] + (1 - mixf) * Ct1[Ncel];       { .. end cell }
      for j := 1 to Ncel do Ct1[j] := Ct2[j];                      { .. copy back into Ct1-cells }
    end;

  writeln(' Depth below seabottom (m); Cl (mol/L) after', Tott/(31536e3):8:2,' years');
  writeln('             0     ;',Ct1[0]:8:4, Nmix:5, ' timesteps');
  for i:= 1 {-3} to Ncel do
    begin    { for 90 cells output }
      i := i {+ 8; if i > Ncel then i := Ncel };
      writeln('          ', Delx * (i-0.5) / 100:9:3, Ct1[i]:10:6);
    end;
end.