Program Columnperc;  { Tritium/Helium dating with variable diffusion coeff's}
Const  Nel = 1;    Ncelmax = 200;
Var Numdisp, bc0_He                         : boolean;
    i, j, cnt                               : integer;
    Ncel, Nshift, Ish, Nmix, Mix_count      : integer;
    Veloc, Tott, Delt, Totx, P              : real;
    Mixf, D, Disp, R, Cort, dc              : real;
    Delx, Age                               : array[1..Ncelmax] of real;
    Cin                                     : array[0..Nel] of real;
    Ct1, Ct2, Mix                        : array[0..Nel, 1..Ncelmax] of real;
    file1, file2                            : text;

procedure dspvty;  			{ calculate mixing factor and dispersive timesteps... }
  begin
    Cort := 1.0; {+ 2 / Ncel; } 	{ Cort corrects for lack of end-cell mixing... }
					{ find max of Mixf and hence Nmix ...}
    Mixf := (D * R * Delt + Disp * Delx[Ncel]) / (Delx[Ncel] * Delx[Ncel]) * Cort;
    if  bc0_He = false then
      Nmix := 1 + trunc(3.0 * Mixf)
      else
      Nmix := 1 + trunc(4.5 * Mixf);
      if Nmix < 2 then Nmix := 2;
					{ define mixing factors for helium,
					  diffuses R times faster than tritium...}
    for j := 1 to Ncel do
      Mix[1, j] := (D * R * Delt + Disp * Delx[j]) / (Delx[j] * Delx[j]) / Nmix * Cort;
					{ define mixing factors for tritium...}
    for j := 1 to Ncel do
      Mix[0, j] := (D * Delt + Disp * Delx[j]) / (Delx[j] * Delx[j]) / Nmix * Cort;
  end;

procedure dffdsp;  			{ diffuse and disperse... }
var Mixup, Mixdn  : real;
  begin
					{  mix inner cells ... }
    for j:= 2 to Ncel-1 do   for i:= 0 to Nel do
					{ average mixing factors of neigboring cells .. }
      begin
        Mixup := (Mix[i, j-1] + Mix[i, j]) / 2;
        Mixdn := (Mix[i, j] + Mix[i, j+1]) / 2;
        if Mixup < 0 then   begin   Mixup := 0;  Numdisp := true;   end;
        if Mixdn < 0 then Mixdn := 0;
        Ct2[i, j] := (1.0 - Mixup - Mixdn) * Ct1[i, j] + Mixup * Ct1[i, j-1] + Mixdn * Ct1[i, j+1];
      end;
					{  mix boundary cells ... }
    for i:= 0 to Nel do
      begin
	if (( bc0_He = false) or (i = 0)) then
	  begin                         { flux bc for tritium ... }
	    Mixdn := (Mix[i, 1] + Mix[i, 2]) / 2;
	    if Mixdn < 0 then Mixdn := 0;
	    Ct2[i, 1] := (1 - Mixdn) * Ct1[i, 1] + Mixdn * Ct1[i, 2];
	  end else
	  if i = 1 then
	    begin                         { constant bc = 0 for helium in cell 1...}
	      Mixdn := (Mix[i, 1] + Mix[i, 2]) / 2;
	      if Mixdn < 0 then Mixdn := 0;
	      Mixup := Mix[i, 1];
	      Ct2[i, 1] := (1 - 2 * Mixup - Mixdn) * Ct1[i, 1] + Mixdn * Ct1[i, 2] + 2 * Cin[i];
	    end;
					{ last cell, closed bc... }
        Mixup := (Mix[i, Ncel-1] + Mix[i, Ncel]) / 2;
        if Mixup < 0 then Mixup := 0;
        Ct2[i, Ncel] := (1 - Mixup) * Ct1[i, Ncel] + Mixup * Ct1[i, Ncel-1];
      end;
    for j:=1 to Ncel do   for i:=0 to Nel do   Ct1[i, j]:=Ct2[i, j];
  end;

procedure flush;                        { advect... }
  begin
    for i:= 0 to Nel do
      begin
	for j := Ncel downto 2 do
          begin
	    Ct1[i, j] := Ct1[i, j-1];   { ... transport }
            Ct2[i, j] := Ct1[i, j];
          end;
         { also the first cell... }
	Ct1[i, 1] := Cin[i];
        Ct2[i, 1] := Ct1[i, 1];
      end;
  end;

procedure decay(j : integer; t1 : real);
  begin
    dc := exp(-0.0558 * (t1 / (1 + Nmix)));
    Ct1[1, j] := Ct1[1, j] + Ct1[0, j] * (1 - dc);
    Ct1[0, j] := Ct1[0, j] * dc;
    Ct2[0, j] := Ct1[0, j];
    Ct2[1, j] := Ct1[1, j];
  end;


{Main}
begin
					{ define column data...}
  Totx := 35.0 { m };  P := 1.0 { m/yr};
  Ncel := 44;				{ change to 100 for Figs 11.9 & 11.10 }
					{ check program with zero diffusion,
					  change to 1.0e-9 in the real world ...}
  D := 0.0e-9 * (3600 * 24.0 * 365)  {m2/yr};
  R := 1;				{ Helium diffuses R = 4 times faster than tritium }
  Disp := 0.0 {m}; 			{ Disp could be 6 m }
  bc0_He := false;			{ set constant concentration boundary true for helium }
  Delt := 1.; { yr }
  Nshift := round(36 / Delt);           { from 1950 ... 1985}
  Delx[1] := Totx * (1 - exp(- Delt * P / Totx));
  Veloc := Delx[1] / Delt {m/yr};
  dc := 0.0;
  for j := 1 to Ncel do
    begin
      Delx[j] := Totx * (1 - exp(- Delt * P * j / Totx)) - dc;
      dc := dc + Delx[j];
    end;
  writeln('Tritium/Helium in groundwater, from 1950 - ', (1949 + round(Nshift * Delt)));
  writeln('  D_e = ', D:5:4, ' m2/yr. ', chr(224), ' = ', Disp:4:2,
	  ' m. Ncel = ', Ncel, '. d_max = ', dc :4:1,
	  ' m. D_profile = ', Totx : 4:1, ' m');
				       { set initial concentrations ... }
  for j := 1 to Ncel do
    begin
      Age[j] := (j - 0.5) * Delt;
      dc := exp(-0.0558 * Age[j]);
      Ct1[0, j] := 5.0 * dc;
      Ct1[1, j] := 5.0 * (1.0 - dc);
    end;
  Cin[0] := 1.;   Cin[1] := 0.0 {TU};
					{ find mixing factors... }
  dspvty;
					{ now the experiment ...}
  assign(file2, 'tr_rain'); reset(file2);
  readln(file2); readln(file2);
  read(file2, Cin[0]);
  Cin[0] := Cin[0] * exp(-0.0558 * 1);  { decay tritium 1 yr in unsaturated zone }
  cnt := 1;
  writeln('  Tritium input from file `tr_rain`. At d = 0 m, c_He = 0...?  ', bc0_He);
  write('  Calculations ---> ');
  for Ish := 1 to Nshift do
    begin
      write(round(1949 + Ish * Delt), '.');
					{read tritium in rain...}
      if round(Ish * Delt) > cnt then
	begin
	  read(file2, Cin[0]);		{ decay 1 year in unsaturated zone...}
          Cin[0] := Cin[0] * exp(-0.0558 * 1);
	  cnt := cnt + 1;
	end;
					{ disperse and decay... }
	for Mix_count := 1 to trunc(Nmix / 2) do
	  begin
	    write(' Mix ', Mix_count:4, ' of ', Nmix:4, '.');
	    dffdsp;
	    for j := 1 to Ncel do decay(j, Delt);
	    for j := 1 to 18 do write(chr(8));
	  end;
					{ advect and decay... }
      flush;
      decay(1, Delt / 2);
      for j := 2 to Ncel do decay(j, Delt);
					{ disperse and decay... }
      for Mix_count := 1 + trunc(Nmix / 2) to Nmix do
	begin
	  write(' Mix ', Mix_count:4, ' of ', Nmix:4, '.');
	  dffdsp;
	  for j := 1 to Ncel do decay(j, Delt);
	  for j := 1 to 18 do write(chr(8));
	end;
      for j := 1 to 5 do write(chr(8));
    end;
    writeln;
  close(file2);

  assign(file1, 'th.txt'); rewrite(file1);
  writeln(file1, 'cell      x/m    age  Tr/He_age/yr      Tr       He');
  Totx := 0.0;
  for j := 1 to Ncel do
    begin
      Age[j] := 0.0;
      if Ct1[0, j] > 0 then Age[j] := ln( Ct1[0, j] / (Ct1[0, j] + Ct1[1, j])) / -0.0558;
      write(file1, j : 6, (Totx + Delx[j] / 2) : 9:3, (j - 0.5) * Delt : 6:2, Age[j] : 10:3, '  ',
		   Ct1[0, j]  : 12:4, Ct1[1, j]  : 12:4);
      writeln(file1);
      Totx := Totx + Delx[j];
    end;
  close(file1);
  writeln('cell    x/m    age  Tr/He_age/yr     Tr          He   (v_H2O = ', Veloc : 5:3, ' m/yr)');
  Totx := 0;
  for j := 1 to Ncel do
    begin
      writeln( j : 4, (Totx + Delx[j] / 2) : 9:3, (j - 0.5) * Delt : 6:2, Age[j]: 9:3, '  ',
	       Ct1[0, j]  : 12:4, Ct1[1, j]  : 12:4);
      Totx := Totx + Delx[j];
    end;
end. { Columnperc}