Program Columnperc;
Const  Nel=1;    Ncelmax=80;
Var    Numdisp                                    : boolean;
       i, j                                       : integer;
       Ncel, Nshift, Ish, Nmix, Mi                : integer;
       Porv, Veloc, Tott, Delt, Totx, Delx, Cort  : real;
       Mixf, D, Disp, R, Rmax                     : real;
       Cin, Kf, nf                                : array[0..Nel] of real;
       Ct1, Ct2, Ctold, Q, Mix                    : array[0..Nel, 1..Ncelmax] of real;
       Pv                                         : array[1..200] of real;
       Cout                                       : array[0..Nel, 1..200] of real;
       file1                                      : text;

procedure dspvty;  { calculate mixing factor and dispersive timesteps... }
  begin
    Cort := 1 + 2 / Ncel;  { Cort corrects for lack of end-cell mixing... }
    Mixf := (D + Disp * Veloc) * Delt / (Delx * Delx) * Cort;
    Nmix := 1 + trunc(3.0 * Mixf);
    Mixf := Mixf / Nmix;
  { correct mixing factors for conservative element...}
    for j := 1 to Ncel do   Mix[0, j] := Mixf - (1 - Veloc * Delt / Delx) / 2 / Nmix * Cort;
  end;

Procedure dspcor;
Var Cdif : real;
  begin  {  Retardation is calculated as R ... }
    for j := 1 to Ncel do   for i := 1 to Nel do
      begin
        Cdif := Ct1[i,j] - Ctold[i,j];
        if abs(Cdif) > 1e-9 then R := (Ct2[i,j] - Ctold[i,j]) / Cdif;
        if R < 1 then R := 1.0;
        if Rmax < R then Rmax := R;
        Mix[i,j] := Mixf - (1 - Veloc * Delt / R / Delx) / 2 / Nmix * Cort;
      end;
  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
        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];
        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... }
var x_in_cell  : real;
  begin
    x_in_cell := Veloc * Delt / Delx;
    for i:= 0 to Nel do
      begin
        for j:=Ncel downto 2 do
          begin
            Ctold[i, j] := Ct1[i, j];   { ... keep old conc's for dspcor }
            Ct1[i, j] := x_in_cell * Ct1[i, j-1] + (1 - x_in_cell) * Ct1[i, j];   { ... transport }
            Ct2[i, j] := Ct1[i, j];
          end;
         { also the first cell... }
        Ctold[i, 1] := Ct1[i, 1];
        Ct1[i, 1] := x_in_cell * Cin[i] + (1 - x_in_cell) * Ct1[i, 1];
        Ct2[i, 1] := Ct1[i, 1];
     end;
  end;

function frndl (Kf, nf, C : real) : real;
  begin
    if C > 1e-10 then frndl := Kf * exp(nf * ln(C)) else frndl := 0;
  end;

procedure distri;
{ distribution of Q and C according to Freundlich-isotherm .. }
Var C1, C2, T, Fc1, Fc2 : real;
    ic, ie, it : integer;
begin
 for ic := 1 to Ncel do  for ie := 1 to Nel do
   begin
     T := Ct1[ie, ic] + Q[ie, ic];
     { .. T(otal) is redistributed over Ct1 and Q with Newton's method ... }
     C2 := T / (Kf[ie] + 1);
     it := 0;
     repeat
       C1 := C2;
       Fc1 := frndl(Kf[ie], nf[ie], C1) + C1 - T;
       C2 := C1 + 1e-9;
       Fc2 := frndl(Kf[ie], nf[ie], C2) + C2 - T;
       if abs(Fc1 - Fc2) > 0 then  C2 := C1 - Fc1 / ((Fc2 - Fc1) / 1e-9);
       if C2 < 1e-9 then C2 := 0.0;
       if C2 > T then C2 := T;
       it := it + 1;
       if it > 20 then begin
         writeln('FREUNDLICH eqn., iteration > 20'); halt; end;
     until (abs(Fc1) <= 1e-6)  or  ((C1 + C2) = 0);
     Ct1[ie, ic] := C2;  Q[ie, ic] := (T - C2);
   end;
end;

procedure fill;
begin
  Delx := Totx / Ncel;
  for i := 0 to Nel do   for j := 1 to Ncel do
    begin
      Ct1[i,j]:=0; Q[i,j]:=0; { .. or, if Ct1>0, then Q=frndl(.,.,.)}
    end;
   Cin[0]:=1;   Cin[1]:=1 {ug/l};
end;

{Main}
begin
  D := 0.0 { cm2/s };   Disp := 0.5 { cm };   Totx := 5.0 { cm };
  Porv := pi * (2.5 * 2.5) * Totx * 0.3 { mL };
  Veloc := 10 / Porv / 3600 * Totx { cm/s };
  Tott := Totx / Veloc * 5; { Time for injection of 5 PV, s};
  Kf[1] := 1; nf[1] := 0.5;
  { initial run for estimating minimum no. of cells ...}
  Ncel := 3;   Delt := Totx / Veloc / 3 { s };   fill;
  Rmax := 1;   dspvty;   Numdisp := false;
  for Ish := 1 to 3 do begin   flush;   distri;   dspcor;   dffdsp;   distri;   end;
  if Numdisp = true then  Ncel := 2 + trunc(Totx / Disp * (1 - 1 / Rmax) / 2);

  { Now start realistic simulation ...}
  fill;   Rmax := 1;   Numdisp := false;
  Nshift := 1 + trunc(Tott * Veloc / Delx);
  Delt := Tott / Nshift {s};
  dspvty;
  for Ish:=1 to Nshift do
    begin
      flush;  { .. new conc's in Ct1 and Ct2; old conc's in Ctold }
      distri; { .. Ct1 has new conc's, after distri }
      dspcor; { .. difference amongst Ct1, Ct2, Ctold gives num. dispersion }
           {    mix, to model remaining dispersion ... }
      for Mi := 1 to Nmix do   begin   dffdsp;   distri;   end;
      Pv[Ish] := (Ish * Delt * Veloc / Delx + 0.5) / Ncel;
      for i:=0 to Nel do Cout[i, Ish]:=Ct1[i, Ncel];
    end;
  assign(file1,'ex10_5'); rewrite(file1);
  for j := 1 to Nshift do
    begin
      write(file1, Pv[j]:6:3,'  ');
      for i := 0 to Nel do write(file1, Cout[i, j]:8:3); writeln(file1);
    end;
  close(file1);
end.