unit erg_color;

interface

procedure SplitColoredElements(tracerBased: Boolean; AutoBurialFluxes: Boolean);

implementation

uses erg_base, erg_types, SysUtils, classes, Dialogs, math;

type tRealArray = Array of Real;
type tRealArray2d = Array of TRealArray;

var newCElements: Array of TErgCElement;

procedure SplitCElements();
// this procedure creates the following tracers:
//   total_red_N, total_red_N_at_bottom
//   aged_red_N, aged_red_N_at_bottom
// and the following processes:
//   aging_red_N, aging_red_N_at_bottom
var i, r, rt, rtr: Integer;
begin
  for i:=0 to length(cElements)-1 do
  begin
    if cElements[i].myIsTracer=1 then
    begin
      //adding 3d tracer
      setLength(tracers,length(tracers)+1);
      initErgTracer(tracers[length(tracers)-1]);
      tracers[length(tracers)-1].name:='total_'+cElements[i].color+'_'+cElements[i].element;
      tracers[length(tracers)-1].description:=tracers[length(tracers)-1].name;
      tracers[length(tracers)-1].isPositive:=0;  //prevent useless mixing, as tracer is calculated every time
      tracers[length(tracers)-1].atmosDep:=0;
      tracers[length(tracers)-1].riverDep:=0;
      tracers[length(tracers)-1].initValue:='0.0';
      tracers[length(tracers)-1].useInitValue:=1;

      //adding bottom tracer
      setLength(tracers,length(tracers)+1);
      initErgTracer(tracers[length(tracers)-1]);
      tracers[length(tracers)-1].name:='total_'+cElements[i].color+'_'+cElements[i].element+'_at_bottom';
      tracers[length(tracers)-1].description:=tracers[length(tracers)-1].name;
      tracers[length(tracers)-1].vertLoc:=1;
      tracers[length(tracers)-1].initValue:='0.0';
      tracers[length(tracers)-1].useInitValue:=1;
      tracers[length(tracers)-1].isPositive:=0;

      //adding top tracer
{      setLength(tracers,length(tracers)+1);
      initErgTracer(tracers[length(tracers)-1]);
      tracers[length(tracers)-1].name:='total_'+cElements[i].color+'_'+cElements[i].element+'_at_top';
      tracers[length(tracers)-1].description:=tracers[length(tracers)-1].name;
      tracers[length(tracers)-1].index:='idx_flat_'+tracers[length(tracers)-1].name;
      tracers[length(tracers)-1].isFlat:=2;
      tracers[length(tracers)-1].rangeMin:='-1.0E+20';
      tracers[length(tracers)-1].rangeMax:='1.0E+20';
      tracers[length(tracers)-1].initValue:='0.0';
      tracers[length(tracers)-1].useInitValue:=1;
      tracers[length(tracers)-1].isPositive:=0;
}    end;

    if cElements[i].myIsAging=1 then
    begin
      //adding 3d tracer
      setLength(tracers,length(tracers)+1);
      initErgTracer(tracers[length(tracers)-1]);
      tracers[length(tracers)-1].name:='aged_'+cElements[i].color+'_'+cElements[i].element;
      tracers[length(tracers)-1].description:=tracers[length(tracers)-1].name;
      tracers[length(tracers)-1].isPositive:=0;  //prevent useless mixing, as tracer is calculated every time
      tracers[length(tracers)-1].atmosDep:=0;
      tracers[length(tracers)-1].riverDep:=0;
      tracers[length(tracers)-1].initValue:='0.0';
      tracers[length(tracers)-1].useInitValue:=1;

      //aging process for 3d tracer
      setLength(processes,length(processes)+1);
      InitErgProcess(processes[length(processes)-1]);
      processes[length(processes)-1].name:='aging_'+cElements[i].color+'_'+cElements[i].element;
      processes[length(processes)-1].description:=processes[length(processes)-1].name;
      processes[length(processes)-1].turnover:='total_'+cElements[i].color+'_'+cElements[i].element;
      SetLength(processes[length(processes)-1].output,1);
      processes[length(processes)-1].output[0].tracer:='aged_'+cElements[i].color+'_'+cElements[i].element;
      processes[length(processes)-1].output[0].amount:='1.0';

      //adding bottom tracer
      setLength(tracers,length(tracers)+1);
      initErgTracer(tracers[length(tracers)-1]);
      tracers[length(tracers)-1].name:='aged_'+cElements[i].color+'_'+cElements[i].element+'_at_bottom';
      tracers[length(tracers)-1].description:=tracers[length(tracers)-1].name;
      tracers[length(tracers)-1].vertLoc:=1;
      tracers[length(tracers)-1].initValue:='0.0';
      tracers[length(tracers)-1].useInitValue:=1;
      tracers[length(tracers)-1].isPositive:=0;

      //aging process for bottom tracer
      setLength(processes,length(processes)+1);
      InitErgProcess(processes[length(processes)-1]);
      processes[length(processes)-1].name:='aging_'+cElements[i].color+'_'+cElements[i].element+'_at_bottom';
      processes[length(processes)-1].description:=processes[length(processes)-1].name;
      processes[length(processes)-1].turnover:='total_'+cElements[i].color+'_'+cElements[i].element+'_at_bottom';
      processes[length(processes)-1].vertLoc:=1;
      SetLength(processes[length(processes)-1].output,1);
      processes[length(processes)-1].output[0].tracer:='aged_'+cElements[i].color+'_'+cElements[i].element+'_at_bottom';
      processes[length(processes)-1].output[0].amount:='1.0';

      GenerateIndexes;
    end;
  end;
end;

procedure SplitElements();
// adds colored elements to the element list
var i: Integer;
begin
  setLength(newCElements,length(cElements));
  for i:=0 to length(cElements)-1 do
  begin
    InitErgCElement(newCElements[i]);
    copyErgCElement(cElements[i],newCElements[i]);
    setLength(elements,length(elements)+1);
    InitErgElement(elements[length(elements)-1]);
    elements[cElements[i].myElementNum].hasColoredCopies:=1;
    CopyErgElement(elements[cElements[i].myElementNum], elements[length(elements)-1]);
    elements[length(elements)-1].name:=cElements[i].color+'_'+elements[length(elements)-1].name;
    elements[length(elements)-1].description:=cElements[i].description;
    elements[length(elements)-1].isColored:=1;
    newCElements[i].element:=elements[length(elements)-1].name;
    newCElements[i].myElementNum:=length(elements)-1;
  end;
end;

procedure SplitTracers(tracerBased: Boolean);
// adds additional tracers like nn_with_red_N and sets their properties
var i, j: Integer;
    t, t_old_max: Integer;
    c: Integer;
    r, rt, rtr: Integer;
    myk,k: Integer;
    gefunden: Boolean;
    isContained: Boolean;
    tracerNameList: TStringList;
    tracerInList: Boolean;
    s: String;
begin
  DecimalSeparator:='.';
  t_old_max:=length(tracers)-1;
  tracerNameList := tStringList.Create;
  for i:=0 to length(cElements)-1 do
  begin
    tracerNameList.Clear;
    if trim(lowercase(CElements[i].duplicatedTracers))<>'all' then
    begin
      s:=trim(lowercase(CElements[i].duplicatedTracers));
      while length(s)>0 do
        tracerNameList.Add(SemiItem(s));
    end;
    for t:=0 to t_old_max do
    begin
      if tracerNameList.Count=0 then
        tracerInList:=true
      else
      begin
        tracerInList:=false;
        for j:=0 to tracerNameList.Count-1 do
          if trim(lowercase(tracers[t].name))=trim(lowercase(tracerNameList[j])) then
            tracerInList:=true;
      end;
      if (tracers[t].isActive>0) and (tracerInList) then
      begin
        isContained:=false;
        for c:=0 to length(tracers[t].contents)-1 do
          if tracers[t].contents[c].myElementNum=cElements[i].myElementNum then
            isContained:=true;
        if IsContained then  //duplicate tracer
        begin
          setLength(tracers,length(tracers)+1);
          initErgTracer(tracers[length(tracers)-1]);
          tracers[length(tracers)-1].comment:=tracers[t].comment;
          tracers[length(tracers)-1].vertSpeed:=tracers[t].vertSpeed;
          tracers[length(tracers)-1].vertDiff:=tracers[t].vertDiff;
          tracers[length(tracers)-1].solubility:=tracers[t].solubility;
          tracers[length(tracers)-1].schmidtNumber:=tracers[t].schmidtNumber;
          tracers[length(tracers)-1].gasName:=tracers[t].gasName;
          if (tracers[t].solubility<>'0') and (tracers[t].gasName='') then //the name of the gas is - per default - the name of the tracer
            tracers[length(tracers)-1].gasName:=tracers[t].name; //then the child tracer has the same gas as the parent
          tracers[length(tracers)-1].vertLoc:=tracers[t].vertLoc;
          tracers[length(tracers)-1].isInPorewater:=tracers[t].isInPorewater;
          tracers[length(tracers)-1].isPositive:=tracers[t].isPositive;
          SetLength(tracers[length(tracers)-1].contents,1);
          tracers[length(tracers)-1].childOf:=tracers[t].name;
          tracers[length(tracers)-1].myChildOf:=t;
          tracers[length(tracers)-1].opacity:='0';
          tracers[length(tracers)-1].atmosDep:=tracers[t].atmosDep*cElements[i].atmosDep;
          tracers[length(tracers)-1].riverDep:=tracers[t].riverDep*cElements[i].riverDep;
          tracers[length(tracers)-1].initValue:='0.0';
          tracers[length(tracers)-1].useInitValue:=1;
          tracers[length(tracers)-1].isPositive:=0;
          tracers[length(tracers)-1].tracerAbove:=tracers[t].tracerAbove;
          tracers[length(tracers)-1].tracerBelow:=tracers[t].tracerBelow;
          tracers[length(tracers)-1].longname:=cElements[i].longname_prefix+tracers[t].longname+cElements[i].longname_suffix;
          tracers[length(tracers)-1].outputUnit:=tracers[t].outputUnit;
          tracers[length(tracers)-1].stdname_prefix:=tracers[t].stdname_prefix;
          tracers[length(tracers)-1].stdname_suffix:=tracers[t].stdname_suffix;
          tracers[length(tracers)-1].myTracerAbove:=-1;
          tracers[length(tracers)-1].myTracerBelow:=-1;
          tracers[t].numChildren:=tracers[t].numChildren+1;

          if tracerBased then
          begin
            tracers[length(tracers)-1].name:=tracers[t].name+'_with_'+cElements[i].color+'_'+celements[i].element;
            if cElements[i].longname_suffix='' then
               tracers[length(tracers)-1].longname:=tracers[length(tracers)-1].longname+'_with_'+cElements[i].color+'_'+celements[i].element;
            tracers[length(tracers)-1].description:=tracers[t].description+'; containing '+cElements[i].description;
            tracers[length(tracers)-1].initValue:='0.0';
            //   We do this further below!
            if tracers[length(tracers)-1].tracerAbove <> 'none' then
              tracers[length(tracers)-1].tracerAbove:=tracers[t].tracerAbove+'_with_'+cElements[i].color+'_'+celements[i].element;
            if tracers[length(tracers)-1].tracerBelow <> 'none' then
              tracers[length(tracers)-1].tracerBelow:=tracers[t].tracerBelow+'_with_'+cElements[i].color+'_'+celements[i].element;

            for c:=0 to length(tracers[t].contents)-1 do
              if tracers[t].contents[c].myElementNum=cElements[i].myElementNum then
              begin
                tracers[length(tracers)-1].contents[0]:=tracers[t].contents[c];
                tracers[length(tracers)-1].contents[0].element:=newCElements[i].element;
                tracers[length(tracers)-1].contents[0].myElementNum:=newCElements[i].myElementNum;
              end;
          end
          else
          begin
            tracers[length(tracers)-1].name:=cElements[i].color+'_'+celements[i].element+'_in_'+tracers[t].name;
            tracers[length(tracers)-1].longname:=cElements[i].color+'_'+celements[i].element+'_in_'+tracers[t].longname;
            tracers[length(tracers)-1].description:=cElements[i].description+' contained in '+tracers[t].description;
            if tracers[length(tracers)-1].tracerAbove <> 'none' then
              tracers[length(tracers)-1].tracerAbove:=cElements[i].color+'_'+celements[i].element+'_in_'+tracers[t].tracerAbove;
            if tracers[length(tracers)-1].tracerBelow <> 'none' then
              tracers[length(tracers)-1].tracerBelow:=cElements[i].color+'_'+celements[i].element+'_in_'+tracers[t].tracerBelow;

            for c:=0 to length(tracers[t].contents)-1 do
              if tracers[t].contents[c].myElementNum=cElements[i].myElementNum then
              begin
                tracers[length(tracers)-1].contents[0]:=tracers[t].contents[c];
                tracers[length(tracers)-1].contents[0].element:=newCElements[i].element;
                tracers[length(tracers)-1].contents[0].myElementNum:=newCElements[i].myElementNum;
                tracers[length(tracers)-1].contents[0].amount:='1';
                tracers[length(tracers)-1].contents[0].myAmount:=1;
                tracers[length(tracers)-1].initValue:='0.0';
              end;
          end;
        end;
      end;
    end;
  end;
  tracerNameList.Free;

  // update some reference indices
  for t:=0 to length(tracers)-1 do
  begin
    if ( (tracers[t].myTracerAbove=-1) and (tracers[t].tracerAbove<>'none') ) then
    begin
      for i:=0 to length(tracers)-1 do
        if tracers[t].tracerAbove=tracers[i].name then
          tracers[t].myTracerAbove:=i;
      if tracers[t].myTracerAbove=-1 then
      // see explanation to tracerBelow below
      begin
        tracers[t].tracerAbove:='none';
        // if tracers[t].childOf<>'none' then
        // begin
        //   tracers[t].tracerAbove:=tracers[tracers[t].myChildOf].tracerAbove;
        //   tracers[t].myTracerAbove:=tracers[tracers[t].myChildOf].myTracerAbove;
        // end
        // else
        // begin
        //   tracers[t].tracerAbove:='none';
        // end;
      end;
    end;
    if ( (tracers[t].myTracerBelow=-1) and (tracers[t].tracerBelow<>'none') ) then
    begin
      for i:=0 to length(tracers)-1 do
        if tracers[t].tracerBelow=tracers[i].name then
          tracers[t].myTracerBelow:=i;
      if tracers[t].myTracerBelow=-1 then
      // there was no tracer, which name is tracerBelow
      //  case 1: tracerBelow==none => myTracerBelow should be -1
      //  case 2: tracer was colored, name==parentTraver_with_COL_ELEMENT and
      //            tracerBelow==parentTracer.tracerBelow_with_COL_ELEMENT
      // Nothing would be to do in case 1. In the past in case 2, tracerBelow
      //  was set to parentTracer.tracerBelow. However, this might lead to
      //  problems when there is a flow into tracerBelow: The parent tracer
      //  includes already this tracer. Therefore, one could accidently let flow
      //  too much of colored tracer into the tracer below.
      begin
        // We overwrite tracerBelow with 'none'
        tracers[t].tracerBelow:='none';
        // Previously
        // if tracers[t].childOf<>'none' then
        // begin
        //   tracers[t].tracerBelow:=tracers[tracers[t].myChildOf].tracerBelow;
        //   tracers[t].myTracerBelow:=tracers[tracers[t].myChildOf].myTracerBelow;
        // end
        // else
        // begin
        //   tracers[t].tracerBelow:='none';
        // end;
      end;
    end;
  end;
end;

procedure PreBalanceProcesses(var A3d, ATop, ABottom: TRealArray2d);
// Determine output minus input of each colored element in all
// original processes, in 3d, top and bottom.
// This is needed in routine PostBalanceProcesses to detect
//    - transfer top <-> 3d, bottom <-> 3d
//    - burial from 3d, top, bottom
// of the colored element
// This is done to auto-generate burial fluxes or, in case of aging tracers,
// to transfer or bury the age as well.
//
// Output variables:
// A3d, ATop, ABottom:
// Array[colored Elements, processes] which
// contains the net creation / destruction of this colored element
// by this process
//
// Repainting processes are considered as burial of the old colored element
// and reset the age to zero. That's why, if a colored element is repainted
// during a process, only the sinks of the old element are taken into
// account.
var
  c,p, i,o,t,e, r,j: Integer;
  bal3dI, balTopI, balBottomI: Real;
  bal3dO, balTopO, balBottomO: Real;
  keepsItsColor: Boolean;
  tracerNameList: TStringList;
  tracerInList: Boolean;
  s: String;
begin
  SetLength(A3d     ,length(cElements),length(processes));
  SetLength(ATop    ,length(cElements),length(processes));
  SetLength(ABottom ,length(cElements),length(processes));
  tracerNameList := TStringList.Create;
  for c:=0 to length(cElements)-1 do
  begin
    tracerNameList.Clear;
    if trim(lowercase(CElements[c].duplicatedTracers))<>'all' then
    begin
      s:=trim(lowercase(CElements[c].duplicatedTracers));
      while length(s)>0 do
        tracerNameList.Add(SemiItem(s));
    end;
    for p:=0 to length(processes)-1 do
    begin
      bal3dI:=0; balTopI:=0; balBottomI:=0;  //initializing balances
      bal3dO:=0; balTopO:=0; balBottomO:=0;
      //considering each element input to process as negative
      for i:=0 to length(processes[p].input)-1 do
      begin
        t:=processes[p].input[i].myTracerNum;
        if tracerNameList.Count=0 then
          tracerInList:=true
        else
        begin
          tracerInList:=false;
          for j:=0 to tracerNameList.Count-1 do
            if trim(lowercase(tracers[t].name))=trim(lowercase(tracerNameList[j])) then
              tracerInList:=true;
        end;
        if (tracers[t].isActive>0) and tracerInList then
        begin
          for e:=0 to length(tracers[t].contents)-1 do
            if tracers[t].contents[e].myElementNum = cElements[c].myElementNum then
              if tracers[t].vertLoc = 0 then
                bal3dI:=bal3dI-processes[p].input[i].myAmount*tracers[t].contents[e].myAmount
              else if tracers[t].vertLoc = 1 then
                balBottomI:=balBottomI-processes[p].input[i].myAmount*tracers[t].contents[e].myAmount
              else if tracers[t].vertLoc = 2 then
                balTopI:=balTopI-processes[p].input[i].myAmount*tracers[t].contents[e].myAmount;
        end;
      end;
      //considering each element output from process as positive
      for o:=0 to length(processes[p].output)-1 do
      begin
        t:=processes[p].output[o].myTracerNum;
        if tracerNameList.Count=0 then
          tracerInList:=true
        else
        begin
          tracerInList:=false;
          for j:=0 to tracerNameList.Count-1 do
            if trim(lowercase(tracers[t].name))=trim(lowercase(tracerNameList[j])) then
              tracerInList:=true;
        end;
        if (tracers[t].isActive>0) and tracerInList then
        begin
          for e:=0 to length(tracers[t].contents)-1 do
            if tracers[t].contents[e].myElementNum = cElements[c].myElementNum then
              if tracers[t].vertLoc = 0 then
                bal3dO:=bal3dO+processes[p].output[o].myAmount*tracers[t].contents[e].myAmount
              else if tracers[t].vertLoc = 1 then
                balBottomO:=balBottomO+processes[p].output[o].myAmount*tracers[t].contents[e].myAmount
              else if tracers[t].vertLoc = 2 then
                balTopO:=balTopO+processes[p].output[o].myAmount*tracers[t].contents[e].myAmount;
        end;
      end;
      //searching for repaints that mean a "burial" of the existing colored element
      for r:=0 to length(processes[p].repaint)-1 do
        if processes[p].repaint[r].myElementNum = cElements[c].myElementNum then
          if ((lowercase(processes[p].repaint[r].oldColor)
               = lowercase(cElements[c].color)) or
              (lowercase(processes[p].repaint[r].oldColor)
               = 'all')) and
              (lowercase(processes[p].repaint[r].newColor)
               <> lowercase(cElements[c].color)) then
          begin
            bal3dO:=0; balTopO:=0; balBottomO:=0;  //no output of this colored element
          end;
      //storing calculated values
      A3d[c,p]:=bal3dI+bal3dO;
      ATop[c,p]:=balTopI+balTopO;
      ABottom[c,p]:=balBottomI+balBottomO;
    end;
  end;
  TracerNameList.Free;
end;

procedure SplitProcesses;
// This procedure creates additonal processes which have colored tracers as input.
// These additional processes are used to generate the source terms of colored tracers.
// Their sink terms are generated in another way.
//
// Each process that has input or output of an element which is colored will be duplicated.
// The rate of the new process is determined by:
// r' = r * Sum_i(a_i*b_j'*flat*[A_i'/A_i]) / Sum_i(a_i*b_j*flat)
// here, r    is the rate of the basic process
//       a_i  is .input[i].amount
//       A_i  is original tracer concentration
//       b_j  is .contents[j].amount of that tracer, where j is the considered element
//       A_i' is colored tracer concentration
//       b_j' is .contents[j].amount of the colored tracer
//       flat is usually 1.0; for vertLoc=SED or vertLoc=SUR tracers this is cgt_cellheight*cgt_density to get the units correct
//       [...] means max(0.0,min(1.0,...))

const epsilon='0.00000000001';
var ce: Integer;
    p, pp, pmax: Integer;
    i, o, j: Integer;
    c: Integer;
    t, tt: Integer;
    contained: Boolean;
    ElInputs, cElInputs: TStringList;     // list of total inputs ; colored element inputs to the currently duplicated process
    ContentString: String;
    tracerNameList: TStringList;
    tracerInList: Boolean;
    s: String;
begin
  DecimalSeparator:='.';
  ElInputs:=TStringList.create;
  cElInputs:=TStringList.create;
  tracerNameList := TStringList.Create;
  pmax:=length(processes)-1;
  for ce:=0 to length(cElements)-1 do //loop over colored elements
  begin
    tracerNameList.Clear;
    if trim(lowercase(CElements[ce].duplicatedTracers))<>'all' then
    begin
      s:=trim(lowercase(CElements[ce].duplicatedTracers));
      while length(s)>0 do
        tracerNameList.Add(SemiItem(s));
    end;
    for p:=0 to pmax do  //loop over old processes
    begin
      //check if it's output contains this element
      contained:=false;
      for o:=0 to length(processes[p].output)-1 do
        for c:=0 to length(tracers[processes[p].output[o].myTracerNum].contents)-1 do
        begin
          if tracerNameList.Count=0 then
            tracerInList:=true
          else
          begin
            tracerInList:=false;
            for j:=0 to tracerNameList.Count-1 do
              if trim(lowercase(tracers[processes[p].output[o].myTracerNum].name))=trim(lowercase(tracerNameList[j])) then
                tracerInList:=true;
          end;
          if tracers[processes[p].output[o].myTracerNum].contents[c].myElementNum = cElements[ce].myElementNum then
            if (tracers[processes[p].output[o].myTracerNum].isActive>0) and tracerInList then
              contained:=true;
        end;
      //also check if input contains this element
      for o:=0 to length(processes[p].input)-1 do
        for c:=0 to length(tracers[processes[p].input[o].myTracerNum].contents)-1 do
        begin
          if tracerNameList.Count=0 then
            tracerInList:=true
          else
          begin
            tracerInList:=false;
            for j:=0 to tracerNameList.Count-1 do
              if trim(lowercase(tracers[processes[p].input[o].myTracerNum].name))=trim(lowercase(tracerNameList[j])) then
                tracerInList:=true;
          end;
          if tracers[processes[p].input[o].myTracerNum].contents[c].myElementNum = cElements[ce].myElementNum then
            if (tracers[processes[p].input[o].myTracerNum].isActive>0) and tracerInList then
              contained:=true;
        end;
      //if contained, the process will be duplicated;
      if contained then
      begin
        //step 1: generate a "colored process" which shall generate all the colored outputs
        setLength(processes,length(processes)+1);
        pp:=length(processes)-1;
        InitErgProcess(processes[pp]);
        CopyErgProcess(processes[p],processes[pp]);
        processes[pp].processType:=processes[p].processType;
        processes[pp].name:=processes[pp].name+'_'+cElements[ce].color+'_'+cElements[ce].element;
        processes[pp].description:=processes[pp].description+'; sub-process for '+cElements[ce].color+' '+elements[cElements[ce].myElementNum].description;
        setLength(processes[pp].output,length(processes[p].output));
        setLength(processes[pp].input,0); //only in output, colored tracers have fixed ratio
        setLength(processes[pp].repaint,0);
        for i:=0 to length(processes[p].output)-1 do
          processes[pp].output[i]:=processes[p].output[i];
        ElInputs.Clear;
        cElInputs.Clear;
        //determining the rate as (old rate) * (colored element input) / (total element input)
        for i:=0 to length(processes[p].input)-1 do
        begin
          if tracerNameList.Count=0 then
            tracerInList:=true
          else
          begin
            tracerInList:=false;
            for j:=0 to tracerNameList.Count-1 do
              if trim(lowercase(tracers[processes[p].input[i].myTracerNum].name))=trim(lowercase(tracerNameList[j])) then
                tracerInList:=true;
          end;
          if (tracers[processes[p].input[i].myTracerNum].isActive>0) and tracerInList then
            for c:=0 to length(tracers[processes[p].input[i].myTracerNum].contents)-1 do
              if tracers[processes[p].input[i].myTracerNum].contents[c].myElementNum = cElements[ce].myElementNum then
              begin
                //total element input
                contentString:='('+processes[p].input[i].amount+')*('+tracers[processes[p].input[i].myTracerNum].contents[c].amount+')';
{                if (tracers[processes[p].input[i].myTracerNum].vertLoc=0) and (processes[p].vertLoc=1) then
                  contentString:=contentString+'/'+cellheightTimesDensity
                else if (tracers[processes[p].input[i].myTracerNum].vertLoc=0) and (processes[p].vertLoc=2) then
                  contentString:=contentString+'/'+cellheightTimesDensity;}
                ElInputs.Add(contentString);
                //colored element input
                //finding new tracer
                tt:=-1;
                for t:=0 to length(tracers)-1 do
                  if (tracers[t].name=tracers[processes[p].input[i].myTracerNum].name+'_with_'+cElements[ce].color+'_'+cElements[ce].element)
                      or (tracers[t].name=cElements[ce].color+'_'+cElements[ce].element+'_in_'+tracers[processes[p].input[i].myTracerNum].name) then
                    tt:=t;
                contentString:='('+processes[p].input[i].amount+')*('+tracers[tt].contents[0].amount+')';
{                if (tracers[processes[p].input[i].myTracerNum].vertLoc=0) and (processes[p].vertLoc=1) then
                  contentString:=contentString+'/'+cellheightTimesDensity
                else if (tracers[processes[p].input[i].myTracerNum].vertLoc=0) and (processes[p].vertLoc=2) then
                  contentString:=contentString+'/'+cellheightTimesDensity;}
                contentString:=contentString+'*max(0.0,min(1.0,'+tracers[tt].name+'/max('+
                  epsilon+','+tracers[processes[p].input[i].myTracerNum].name+')))';
                cElInputs.Add(contentString);
              end;
        end;
        if cElInputs.Count=0 then // no input of colored elements -> deleting process
        begin
          processes[length(processes)-1].Free;
          setLength(processes,length(processes)-1)
        end
        else
        begin
          ContentString:='('+cElInputs[0];
          for i:=1 to cElInputs.Count-1 do
            ContentString:=ContentString+'+'+cElInputs[i];
          ContentString:=ContentString+') / ('+ElInputs[0];
          for i:=1 to ElInputs.Count-1 do
            ContentString:=ContentString+'+'+ElInputs[i];
          ContentString:=processes[p].name+' * '+ContentString+')';
          processes[pp].turnover:=ContentString;
          for o:=0 to length(processes[pp].output)-1 do
          begin
            tt:=-1;
            for t:=0 to length(tracers)-1 do
              if (tracers[t].name=tracers[processes[p].output[o].myTracerNum].name+'_with_'+cElements[ce].color+'_'+cElements[ce].element)
                  or (tracers[t].name=cElements[ce].color+'_'+cElements[ce].element+'_in_'+tracers[processes[p].output[o].myTracerNum].name) then
                tt:=t;
            if tt=-1 then //this output result is not colored -> prepare for deleting
              processes[pp].output[o].tracer:=''
            else
            begin
              processes[pp].output[o].tracer:=tracers[tt].name;
              processes[pp].output[o].myTracerNum:=tt;
            end;
          end;
          //now delete unused tracers
          for o:=length(processes[pp].output)-1 downto 0 do
          begin
            if processes[pp].output[o].tracer='' then
            begin
              processes[pp].output[o]:=processes[pp].output[length(processes[pp].output)-1];
              setLength(processes[pp].output,length(processes[pp].output)-1);
            end;
          end;
        end;
        //step 2: if process is at bottom, and it has both colored tracer input from sediment and from water column,
        //        we need to know how much of which compartment to change the age tracer in a correct way.
        //That's why for vertLoc=SED and vertLoc=SUR processes and cElements with isAging=1, we define two more
        //"virtual processes", each with zero input and output, whose process rate scaled with the input ratio of
        //water column and sediment/surface colored elements, respectively.
        if (processes[p].vertLoc>0) and (celements[ce].isAging='1') then
        begin
          //first, water column
          setLength(processes,length(processes)+1);
          pp:=length(processes)-1;
          InitErgProcess(processes[pp]);
          CopyErgProcess(processes[p],processes[pp]);
          processes[pp].processType:=processes[p].processType;
          processes[pp].name:=processes[pp].name+'_'+cElements[ce].color+'_'+cElements[ce].element+'_in_water';
          processes[pp].description:=processes[pp].description+'; dummy sub-process for age determination of '+cElements[ce].color+' '+elements[cElements[ce].myElementNum].description;
          setLength(processes[pp].output,0); //dummy process, does nothing, only its rate needs to be determined
          setLength(processes[pp].input,0);
          setLength(processes[pp].repaint,0);
          ElInputs.Clear;
          cElInputs.Clear;
          //determining the rate as (old rate) * (colored element input) / (total element input)
          for i:=0 to length(processes[p].input)-1 do
          begin
            if tracerNameList.Count=0 then
              tracerInList:=true
            else
            begin
              tracerInList:=false;
              for j:=0 to tracerNameList.Count-1 do
                if trim(lowercase(tracers[processes[p].input[i].myTracerNum].name))=trim(lowercase(tracerNameList[j])) then
                  tracerInList:=true;
            end;
            if (tracers[processes[p].input[i].myTracerNum].isActive>0) and tracerInList then
              for c:=0 to length(tracers[processes[p].input[i].myTracerNum].contents)-1 do
                if (tracers[processes[p].input[i].myTracerNum].contents[c].myElementNum = cElements[ce].myElementNum) and (tracers[processes[p].input[i].myTracerNum].vertLoc=0) then //only search for tracers with vertLoc=WAT which contain the correct element
                begin
                  //total element input
                  contentString:='('+processes[p].input[i].amount+')*('+tracers[processes[p].input[i].myTracerNum].contents[c].amount+')';
                  ElInputs.Add(contentString);
                  //colored element input
                  //finding new tracer
                  tt:=-1;
                  for t:=0 to length(tracers)-1 do
                    if (tracers[t].name=tracers[processes[p].input[i].myTracerNum].name+'_with_'+cElements[ce].color+'_'+cElements[ce].element)
                        or (tracers[t].name=cElements[ce].color+'_'+cElements[ce].element+'_in_'+tracers[processes[p].input[i].myTracerNum].name) then
                      tt:=t;
                  contentString:='('+processes[p].input[i].amount+')*('+tracers[tt].contents[0].amount+')';
                  contentString:=contentString+'*max(0.0,min(1.0,'+tracers[tt].name+'/max('+
                    epsilon+','+tracers[processes[p].input[i].myTracerNum].name+')))';
                  cElInputs.Add(contentString);
                end;
          end;
          if cElInputs.Count=0 then // no input of colored elements in water -> deleting process
          begin
            processes[length(processes)-1].Free;
            setLength(processes,length(processes)-1)
          end
          else
          begin
            ContentString:='('+cElInputs[0];
            for i:=1 to cElInputs.Count-1 do
              ContentString:=ContentString+'+'+cElInputs[i];
            ContentString:=ContentString+') / ('+ElInputs[0];
            for i:=1 to ElInputs.Count-1 do
              ContentString:=ContentString+'+'+ElInputs[i];
            ContentString:=processes[p].name+' * '+ContentString+')';
            processes[pp].turnover:=ContentString;
          end;
          //second, sediment or surface tracers
          setLength(processes,length(processes)+1);
          pp:=length(processes)-1;
          InitErgProcess(processes[pp]);
          CopyErgProcess(processes[p],processes[pp]);
          processes[pp].processType:=processes[p].processType;
          if processes[p].vertLoc=1 then
            processes[pp].name:=processes[pp].name+'_'+cElements[ce].color+'_'+cElements[ce].element+'_at_bottom'
          else
            processes[pp].name:=processes[pp].name+'_'+cElements[ce].color+'_'+cElements[ce].element+'_at_top';
          processes[pp].description:=processes[pp].description+'; dummy sub-process for age determination of '+cElements[ce].color+' '+elements[cElements[ce].myElementNum].description;
          setLength(processes[pp].output,0); //dummy process, does nothing, only its rate needs to be determined
          setLength(processes[pp].input,0);
          setLength(processes[pp].repaint,0);
          ElInputs.Clear;
          cElInputs.Clear;
          //determining the rate as (old rate) * (colored element input) / (total element input)
          for i:=0 to length(processes[p].input)-1 do
          begin
            if tracerNameList.Count=0 then
              tracerInList:=true
            else
            begin
              tracerInList:=false;
              for j:=0 to tracerNameList.Count-1 do
                if trim(lowercase(tracers[processes[p].input[i].myTracerNum].name))=trim(lowercase(tracerNameList[j])) then
                  tracerInList:=true;
            end;
            if (tracers[processes[p].input[i].myTracerNum].isActive>0) and TracerInList then
              for c:=0 to length(tracers[processes[p].input[i].myTracerNum].contents)-1 do
                if (tracers[processes[p].input[i].myTracerNum].contents[c].myElementNum = cElements[ce].myElementNum) and (tracers[processes[p].input[i].myTracerNum].vertLoc=processes[p].vertLoc) then //only search for tracers with vertLoc=SED or SUR which contain the correct element
                begin
                  //total element input
                  contentString:='('+processes[p].input[i].amount+')*('+tracers[processes[p].input[i].myTracerNum].contents[c].amount+')';
                  ElInputs.Add(contentString);
                  //colored element input
                  //finding new tracer
                  tt:=-1;
                  for t:=0 to length(tracers)-1 do
                    if (tracers[t].name=tracers[processes[p].input[i].myTracerNum].name+'_with_'+cElements[ce].color+'_'+cElements[ce].element)
                        or (tracers[t].name=cElements[ce].color+'_'+cElements[ce].element+'_in_'+tracers[processes[p].input[i].myTracerNum].name) then
                      tt:=t;
                  contentString:='('+processes[p].input[i].amount+')*('+tracers[tt].contents[0].amount+')';
                  contentString:=contentString+'*max(0.0,min(1.0,'+tracers[tt].name+'/max('+
                    epsilon+','+tracers[processes[p].input[i].myTracerNum].name+')))';
                  cElInputs.Add(contentString);
                end;
          end;
          if cElInputs.Count=0 then // no input of colored elements in water -> deleting process
          begin
            processes[length(processes)-1].Free;
            setLength(processes,length(processes)-1)
          end
          else
          begin
            ContentString:='('+cElInputs[0];
            for i:=1 to cElInputs.Count-1 do
              ContentString:=ContentString+'+'+cElInputs[i];
            ContentString:=ContentString+') / ('+ElInputs[0];
            for i:=1 to ElInputs.Count-1 do
              ContentString:=ContentString+'+'+ElInputs[i];
            ContentString:=processes[p].name+' * '+ContentString+')';
            processes[pp].turnover:=ContentString;
          end;
        end;
      end;
    end;
  end;
  tracerNameList.Free;
end;

procedure DeleteProcess(pp: Integer);
//deletes one process and maintains the order of the others
var p: Integer;
begin
  setLength(processes[pp].input,0);
  setLength(processes[pp].output,0);
  setLength(processes[pp].repaint,0);
  for p:=pp to length(processes)-2 do
    CopyErgProcess(processes[p+1],processes[p]);
  processes[length(processes)-1].Free;
  setLength(processes,length(processes)-1);
end;

procedure ApplyRepaint(pp, r, e: Integer);
// For a specific original process and a specific element, the coloring process is done.
// This means that the colored subprocess may no longer have the same color
// as output which it had as input.
// So maybe its output needs to be changed to the correct color, or, if it is
// of color "none", we may a) delete the marked process, if its age is not considered
//                         b) keep it and delete its output, if we need it for age burial
var p, ce, ce2, o, omax, c, t: Integer;
    colorstring: String;
    cAmount, tAmount: Real;
begin
  if lowercase(processes[pp].repaint[r].oldColor)='all' then
  begin
    //if all of this element is painted in the same color, we need no sub-processes.
    //instead, the original process will have colored output in the same amount as total output
    //marking all subprocesses for deleting
    for p:=0 to length(processes)-1 do
      for ce:=0 to length(cElements)-1 do
        if cElements[ce].myElementNum=e then
          if processes[p].name = processes[pp].name+'_'+cElements[ce].color+'_'+cElements[ce].element then
            if cElements[ce].isAging='0' then
              processes[p].name:=''
            else
              setLength(processes[p].output,0);
    //delete marked subprocesses
    for p:=length(processes)-1 downto 0 do
      if processes[p].name='' then
        DeleteProcess(p);
    //check if there is a colored element of the desired new color
    colorstring:='none';
    for ce:=0 to length(cElements)-1 do
      if (lowercase(processes[pp].repaint[r].newColor)=lowercase(cElements[ce].color)) and
         (cElements[ce].myElementNum=e) then
        colorstring:=cElements[ce].color+'_'+cElements[ce].element;
    //if this colored element exists, add colored tracers to the standard process output
    //(if it not exists, the "repaint" may be to color "none" which means we do not need colored output)
    if colorstring <> 'none' then
    begin
      omax:=length(processes[pp].output)-1;
      for o:=0 to omax do
        for t:=0 to length(tracers)-1 do
          if tracers[t].myChildOf = processes[pp].output[o].myTracerNum then
            if pos(colorstring,tracers[t].name) >= 1 then //this tracer has to be inserted
            begin
              //determining relationship of element contents -> tracer- or element-based?
              cAmount:=tracers[t].contents[0].myAmount;
              tAmount:=0.0;
              for c:=0 to length(tracers[processes[pp].output[o].myTracerNum].contents)-1 do
                if tracers[processes[pp].output[o].myTracerNum].contents[c].myElementNum=e then
                  tAmount:=tracers[processes[pp].output[o].myTracerNum].contents[c].myAmount;
              //duplicate output
              setLength(processes[pp].output,length(processes[pp].output)+1);
              processes[pp].output[length(processes[pp].output)-1].tracer:=tracers[t].name;
              processes[pp].output[length(processes[pp].output)-1].myTracerNum:=t;
              if cAmount = tAmount then
              begin
                processes[pp].output[length(processes[pp].output)-1].amount:=processes[pp].output[o].amount;
                processes[pp].output[length(processes[pp].output)-1].myAmount:=processes[pp].output[o].myAmount;
              end
              else
              begin
                processes[pp].output[length(processes[pp].output)-1].amount:=processes[pp].output[o].amount+'*'+FloatToStr(TAmount/CAmount);
                processes[pp].output[length(processes[pp].output)-1].myAmount:=processes[pp].output[o].myAmount*tamount/camount;
              end;
            end;
    end;
  end
  else if lowercase(processes[pp].repaint[r].oldColor) = 'none' then  // we only paint the unmarked part
  begin
    //find new marked element and store it in colorstring
    colorstring:='none';
    for ce:=0 to length(cElements)-1 do
      if (lowercase(processes[pp].repaint[r].newColor)=lowercase(cElements[ce].color)) and
         (cElements[ce].myElementNum=e) then
        colorstring:=cElements[ce].color+'_'+cElements[ce].element;
    if colorstring <> 'none' then
    begin
      //we will first paint all of the output in the new color, and then for each existing color we subtract the output
      //step 1: paint all of the element which passes the process, that is add colored elements to the standard process output
      omax:=length(processes[pp].output)-1;
      for o:=0 to omax do
        for t:=0 to length(tracers)-1 do
          if tracers[t].myChildOf = processes[pp].output[o].myTracerNum then
            if pos(colorstring,tracers[t].name) >= 1 then //this tracer has to be inserted
            begin
              //determining relationship of element contents -> tracer- or element-based?
              cAmount:=tracers[t].contents[0].myAmount;
              tAmount:=0.0;
              for c:=0 to length(tracers[processes[pp].output[o].myTracerNum].contents)-1 do
                if tracers[processes[pp].output[o].myTracerNum].contents[c].myElementNum=e then
                  tAmount:=tracers[processes[pp].output[o].myTracerNum].contents[c].myAmount;
              //duplicate output
              setLength(processes[pp].output,length(processes[pp].output)+1);
              processes[pp].output[length(processes[pp].output)-1].tracer:=tracers[t].name;
              processes[pp].output[length(processes[pp].output)-1].myTracerNum:=t;
              if cAmount = tAmount then
              begin
                processes[pp].output[length(processes[pp].output)-1].amount:=processes[pp].output[o].amount;
                processes[pp].output[length(processes[pp].output)-1].myAmount:=processes[pp].output[o].myAmount;
              end
              else
              begin
                processes[pp].output[length(processes[pp].output)-1].amount:=processes[pp].output[o].amount+'*'+FloatToStr(TAmount/CAmount);
                processes[pp].output[length(processes[pp].output)-1].myAmount:=processes[pp].output[o].myAmount*tamount/camount;
              end;
            end;
      //step 2: all subprocesses for colored tracers get an "input" of the destination color which counteracts their painting with the new color
      for p:=0 to length(processes)-1 do
        for ce:=0 to length(cElements)-1 do
          if cElements[ce].myElementNum=e then
            if processes[p].name = processes[pp].name+'_'+cElements[ce].color+'_'+cElements[ce].element then
            begin
              setLength(processes[p].input,length(processes[p].output));
              for o:=0 to length(processes[p].output)-1 do
              begin
                processes[p].input[o].tracer:=processes[p].output[o].tracer;
                for ce2:=0 to length(celements)-1 do  // input of new color is needed, no matter which color the output may have due to repainting
                  if cElements[ce2].element=cElements[ce].element then
                    processes[p].input[o].tracer:=StringReplace(processes[p].input[o].tracer,cElements[ce2].color+'_'+cElements[ce2].element,processes[pp].repaint[r].newColor+'_'+cElements[ce].element,[rfReplaceAll, rfIgnoreCase]);
                processes[p].input[o].amount:=processes[p].output[o].amount;
                for t:=0 to length(tracers)-1 do
                  if trim(lowercase(processes[p].input[o].tracer))=trim(lowercase(tracers[t].name)) then
                    processes[p].input[o].myTracerNum:=t;
                processes[p].input[o].myAmount:=processes[p].output[o].myAmount;
              end;
            end;
    end;
  end
  else
  begin //oldColor <> 'all'/'none', that is, we have a repaint of a specific color
    for p:=0 to length(processes)-1 do
      for ce:=0 to length(cElements)-1 do
        if (cElements[ce].myElementNum=e) and (lowercase(cElements[ce].color)=lowercase(processes[pp].repaint[r].oldColor))then
        begin
          colorstring:=cElements[ce].color+'_'+cElements[ce].element;
          if processes[p].name = processes[pp].name+'_'+cElements[ce].color+'_'+cElements[ce].element then
          // now we found process p which would initially output tracers of the old color
            if lowercase(processes[pp].repaint[r].newColor)='none' then
            begin
              if cElements[ce].isAging='0' then
                processes[p].name:='' //mark for deleting
              else
                setLength(processes[p].output,0);
            end
            else //the output of process p needs another color
            begin
              for o:=0 to length(processes[p].output)-1 do
                if pos(lowercase(colorstring),lowercase(processes[p].output[o].tracer))>0 then //this is a tracer of old color
                begin
                  //change its name and number
                  processes[p].output[o].tracer:=StringReplace(processes[p].output[o].tracer,colorstring,processes[pp].repaint[r].newColor+'_'+cElements[ce].element,[rfReplaceAll, rfIgnoreCase]);
                  processes[p].output[o].myTracerNum:=-1;
                  for t:=0 to length(tracers)-1 do
                    if lowercase(tracers[t].name)=lowercase(processes[p].output[o].tracer) then
                      processes[p].output[o].myTracerNum:=t;
                  if processes[p].output[o].myTracerNum=-1 then //tracer of new color does not exist -> deleting process
                    if cElements[ce].isAging='0' then
                      processes[p].name:='' //mark for deleting
                    else
                      setLength(processes[p].output,0);
                end;
            end;
        end;
    for p:=length(processes)-1 downto 0 do
      if processes[p].name='' then
        DeleteProcess(p);
  end;
end;

procedure ApplyRepaints;
// loops over all processes and calls ApplyRepaint if they have a repaint statement for an element or for all elements
var p, r, e: Integer;
begin
  for p:=0 to length(processes)-1 do
  begin
    if p<length(processes) then     // some processes may have been deleted during the loop
      for r:=0 to length(processes[p].repaint)-1 do
        for e:=0 to length(elements)-1 do
          if pos('_',elements[e].name)<1 then //old element, not a painted one
            if (processes[p].repaint[r].element=elements[e].name) or (lowercase(processes[p].repaint[r].element)='all') then
              ApplyRepaint(p,r,e);
  end;
end;

procedure PostBalanceProcesses(var A3d, ATop, ABottom: TRealArray2d; AutoBurialFluxes: Boolean);
// This procedure generates processes which transfer the age concentration
// (mean age times concentration) from water to sediment and vice versa.
// Also, it creates "age burial processes" if a colored element is removed
// by an unbalanced process, or by a process which repaints it.
const epsil=1.0e-10;
var
  c, p,pp, pmax,ap, f: Integer;
  TTopTo3d, T3dToTop, TBottomTo3d, T3dToBottom: TStringList;
  amount: Real;
  found: boolean;
  s, s1: String;
  tracerNameList: TStringList;
  tracerInList: Boolean;
begin
  // In the previously called procedure PreBalanceProcesses, we obtained the arrays
  // A3d, ATop and ABottom which calculate the net production / consumption of an
  // element by a process, in water, sediment and surface.
  // Now we first try to link e.g. a production in water and a consumption in
  // sediment together, so we can interpret it as a transfer of element from
  // sediment to water.
  // In this case, the age concentration needs to be transferred as well.
  DecimalSeparator:='.';
  tracerNameList:=TStringList.Create;
  for c:=0 to length(cElements)-1 do
  begin
    TracerNameList.Clear;
    if trim(lowercase(CElements[c].duplicatedTracers))<>'all' then
    begin
      s1:=trim(lowercase(CElements[c].duplicatedTracers));
      while length(s1)>0 do
        tracerNameList.Add(SemiItem(s1));
    end;
    TTopTo3d:=TStringList.Create;
    T3dToTop:=TStringList.Create;
    TBottomTo3d:=TStringList.Create;
    T3dToBottom:=TStringList.Create;
    //collecting rates for transfer of age: top <-> 3d <-> bottom
    //writing to StringLists
    for p:=0 to length(processes)-1 do
     if processes[p].vertLoc<>0 then //transfer can only take place at top or bottom
      if pos('_'+cElements[c].color+'_'+cElements[c].element,processes[p].name)>0 then
        for pp:=0 to p-1 do
          if processes[p].name = processes[pp].name+'_'+cElements[c].color+'_'+cElements[c].element then
          begin  //pp is original process, p is sub-process
            if (A3d[c,pp]>epsil) and (ATop[c,pp]<-epsil) then
            begin
              amount := min(abs(A3d[c,pp]), abs(ATop[c,pp]));
              TTopTo3d.Add(processes[p].name+'*'+FloatToStr(amount));
              A3d[c,pp]:=A3d[c,pp]-amount;
              ATop[c,pp]:=ATop[c,pp]+amount;
            end;
            if (A3d[c,pp]<-epsil) and (ATop[c,pp]>epsil) then
            begin
              amount := min(abs(A3d[c,pp]), abs(ATop[c,pp]));
              T3dToTop.Add(processes[p].name+'*'+FloatToStr(amount));
              A3d[c,pp]:=A3d[c,pp]+amount;
              ATop[c,pp]:=ATop[c,pp]-amount;
            end;
            if (A3d[c,pp]>epsil) and (ABottom[c,pp]<-epsil) then
            begin
              amount := min(abs(A3d[c,pp]), abs(ABottom[c,pp]));
              TBottomTo3d.Add(processes[p].name+'*'+FloatToStr(amount));
              A3d[c,pp]:=A3d[c,pp]-amount;
              ABottom[c,pp]:=ABottom[c,pp]+amount;
            end;
            if (A3d[c,pp]<-epsil) and (ABottom[c,pp]>epsil) then
            begin
              amount := min(abs(A3d[c,pp]), abs(ABottom[c,pp]));
              T3dToBottom.Add(processes[p].name+'*'+FloatToStr(amount));
              A3d[c,pp]:=A3d[c,pp]+amount;
              ABottom[c,pp]:=ABottom[c,pp]-amount;
            end;
          end;
    // The StringLists T3dToBottom, ... contain the transfer of elements between
    // top, water, sediment
    // If they are not zero, we need to create an age concentration transfer process.
    //
    // For conservative processes, the A3d, ABottom and ATop values should be zero now.
    // If processes are not conservative, but consume an aging colored element,
    // we have to include a burial process for the age tracer
    if cElements[c].myIsAging=1 then
    begin
      //generating age-transfer processes if stringList is not empty
      if TTopTo3d.Count>0 then
      begin
        p:=length(processes);
        setLength(processes,p+1);
        InitErgProcess(processes[p]);
        processes[p].name:='t_age_'+cElements[c].color+'_'+cElements[c].element+'_top_3d';
        processes[p].description:='age transfer of '+cElements[c].description+' from top to 3d';
        processes[p].vertLoc:=2;
        setLength(processes[p].input,1);
        processes[p].input[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element+'_at_top';
        processes[p].input[0].amount:='1.0';
        setLength(processes[p].output,1);
        processes[p].output[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element;
        processes[p].output[0].amount:='1.0';
        processes[p].turnover:='max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+'_at_top) / max(total_'+cElements[c].color+'_'+cElements[c].element+'_at_top, '+ModelInfos.ageEpsilon+') * ((';
        for pp:=0 to TTopTo3d.Count-2 do
          processes[p].turnover:=processes[p].turnover+TTopTo3d[pp]+')+(';
        for pp:=TTopTo3d.Count-1 to TTopTo3d.Count-1 do
          processes[p].turnover:=processes[p].turnover+TTopTo3d[pp]+'))';
      end;
      if T3dToTop.Count>0 then
      begin
        p:=length(processes);
        setLength(processes,p+1);
        InitErgProcess(processes[p]);
        processes[p].name:='t_age_'+cElements[c].color+'_'+cElements[c].element+'_3d_top';
        processes[p].description:='age transfer of '+cElements[c].description+' from 3d to top';
        processes[p].vertLoc:=2;
        setLength(processes[p].input,1);
        processes[p].input[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element;
        processes[p].input[0].amount:='1.0';
        setLength(processes[p].output,1);
        processes[p].output[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element+'_at_top';
        processes[p].output[0].amount:='1.0';
        processes[p].turnover:='max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+') / max(total_'+cElements[c].color+'_'+cElements[c].element+', '+ModelInfos.ageEpsilon+') * ((';
        for pp:=0 to T3dToTop.Count-2 do
          processes[p].turnover:=processes[p].turnover+T3dToTop[pp]+')+(';
        for pp:=T3dToTop.Count-1 to T3dToTop.Count-1 do
          processes[p].turnover:=processes[p].turnover+T3dToTop[pp]+'))';
      end;
      if TBottomTo3d.Count>0 then
      begin
        p:=length(processes);
        setLength(processes,p+1);
        InitErgProcess(processes[p]);
        processes[p].name:='t_age_'+cElements[c].color+'_'+cElements[c].element+'_bottom_3d';
        processes[p].description:='age transfer of '+cElements[c].description+' from bottom to 3d';
        processes[p].vertLoc:=1;
        setLength(processes[p].input,1);
        processes[p].input[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom';
        processes[p].input[0].amount:='1.0';
        setLength(processes[p].output,1);
        processes[p].output[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element;
        processes[p].output[0].amount:='1.0';
        processes[p].turnover:='max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom) / max(total_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom, '+ModelInfos.ageEpsilon+') * ((';
        for pp:=0 to TBottomTo3d.Count-2 do
          processes[p].turnover:=processes[p].turnover+TBottomTo3d[pp]+')+(';
        for pp:=TBottomTo3d.Count-1 to TBottomTo3d.Count-1 do
          processes[p].turnover:=processes[p].turnover+TBottomTo3d[pp]+'))';
      end;
      if T3dToBottom.Count>0 then
      begin
        p:=length(processes);
        setLength(processes,p+1);
        InitErgProcess(processes[p]);
        processes[p].name:='t_age_'+cElements[c].color+'_'+cElements[c].element+'_3d_bottom';
        processes[p].description:='age transfer of '+cElements[c].description+' from 3d to bottom';
        processes[p].vertLoc:=1;
        setLength(processes[p].input,1);
        processes[p].input[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element;
        processes[p].input[0].amount:='1.0';
        setLength(processes[p].output,1);
        processes[p].output[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom';
        processes[p].output[0].amount:='1.0';
        processes[p].turnover:='max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+') / max(total_'+cElements[c].color+'_'+cElements[c].element+', '+ModelInfos.ageEpsilon+') * ((';
        for pp:=0 to T3dToBottom.Count-2 do
          processes[p].turnover:=processes[p].turnover+T3dToBottom[pp]+')+(';
        for pp:=T3dToBottom.Count-1 to T3dToBottom.Count-1 do
          processes[p].turnover:=processes[p].turnover+T3dToBottom[pp]+'))';
      end;
    end;
    //generating burial fluxes and age burial processes for colored elements
    pmax:=length(processes)-1;
    for p:=0 to pmax do
      if pos('_'+cElements[c].color+'_'+cElements[c].element,processes[p].name)>0 then
        for pp:=0 to p-1 do
          if processes[p].name = processes[pp].name+'_'+cElements[c].color+'_'+cElements[c].element then
          begin  //pp is original process, p is sub-process
           if (processes[pp].isActive<>0) then
           begin
           found:=false;
           s:=modelinfos.inactiveProcessTypes;
           while s<>'' do
           begin
             s1:=semiItem(s);
             if trim(lowercase(s1))=trim(lowercase(processes[pp].processType)) then found:=true;
           end;
           if found=false then
           begin
            if (A3d[c,pp]<-epsil) then //3d burial
            begin
              if cElements[c].myIsAging=1 then //generate burial process for the age
              begin
                ap:=length(processes);
                setLength(processes,ap+1);
                InitErgProcess(processes[ap]);
                processes[ap].name:='bur_aged_'+cElements[c].color+'_'+cElements[c].element+'_'+processes[pp].name;
                processes[ap].description:='Burial of aged '+cElements[c].description+' by '+processes[pp].description+' in water column';
                processes[ap].vertLoc:=processes[pp].vertLoc;
                if processes[ap].vertLoc=0 then
                  processes[ap].turnover:=FloatToStr(-A3d[c,pp])+' * '+processes[p].name+
                                      ' * max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+
                                      ') / max(total_'+cElements[c].color+'_'+cElements[c].element
                                       +', '+ModelInfos.ageEpsilon+')'
                else if processes[ap].vertLoc<>0 then //sediment or surface process with 3d burial -> must contain 3d input
                  processes[ap].turnover:=FloatToStr(-A3d[c,pp])+' * '+processes[p].name+'_in_water'+
                                      ' * max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+
                                      ') / max(total_'+cElements[c].color+'_'+cElements[c].element
                                       +', '+ModelInfos.ageEpsilon+')';
                SetLength(processes[ap].input,1);
                processes[ap].input[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element;
                processes[ap].input[0].amount:='1.0';
              end;
            end;
            if (ATop[c,pp]<-epsil) then //Top burial
            begin
              if cElements[c].myIsAging=1 then //generate burial process for the age
              begin
                ap:=length(processes);
                setLength(processes,ap+1);
                InitErgProcess(processes[ap]);
                processes[ap].name:='bur_aged_'+cElements[c].color+'_'+cElements[c].element+'_'+processes[pp].name+'_at_top';
                processes[ap].description:='Burial of aged '+cElements[c].description+' by '+processes[pp].description+' at top';
                processes[ap].vertLoc:=processes[pp].vertLoc;
                if processes[ap].vertLoc=0 then
                  processes[ap].turnover:=FloatToStr(-ATop[c,pp])+' * '+processes[p].name+
                                      ' * max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+'_at_top)'+
                                      ' / max(total_'+cElements[c].color+'_'+cElements[c].element+'_at_top'
                                       +', '+ModelInfos.ageEpsilon+')'
                else if processes[ap].vertLoc<>0 then //surface process with surface burial -> must contain surface input
                  processes[ap].turnover:=FloatToStr(-ATop[c,pp])+' * '+processes[p].name+'_at_top'+
                                      ' * max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+'_at_top)'+
                                      ' / max(total_'+cElements[c].color+'_'+cElements[c].element+'_at_top'
                                       +', '+ModelInfos.ageEpsilon+')';
                SetLength(processes[ap].input,1);
                processes[ap].input[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element+'_at_top';
                processes[ap].input[0].amount:='1.0';
              end;
            end;
            if (ABottom[c,pp]<-epsil) then //Bottom burial
            begin
              if cElements[c].myIsAging=1 then //generate burial process for the age
              begin
                ap:=length(processes);
                setLength(processes,ap+1);
                InitErgProcess(processes[ap]);
                processes[ap].name:='bur_aged_'+cElements[c].color+'_'+cElements[c].element+'_'+processes[pp].name+'_at_bottom';
                processes[ap].description:='Burial of aged '+cElements[c].description+' by '+processes[pp].description+' at bottom';
                processes[ap].vertLoc:=processes[pp].vertLoc;
                if processes[ap].vertLoc=0 then
                  processes[ap].turnover:=FloatToStr(-ABottom[c,pp])+' * '+processes[p].name+
                                      ' * max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom)'+
                                      ' / max(total_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom'
                                       +', '+ModelInfos.ageEpsilon+')'
                else if processes[ap].vertLoc<>0 then //surface process with surface burial -> must contain surface input
                  processes[ap].turnover:=FloatToStr(-ABottom[c,pp])+' * '+processes[p].name+'_at_bottom'+
                                      ' * max(0.0,aged_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom)'+
                                      ' / max(total_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom'
                                       +', '+ModelInfos.ageEpsilon+')';
                SetLength(processes[ap].input,1);
                processes[ap].input[0].tracer:='aged_'+cElements[c].color+'_'+cElements[c].element+'_at_bottom';
                processes[ap].input[0].amount:='1.0';
              end;
            end;
           end;
           end;
          end;
  end;
  TracerNameList.Free;
end;

procedure SplitColoredElements(tracerBased: Boolean; AutoBurialFluxes: Boolean);
var A3d, ATop, ABottom: TRealArray2d;
begin
  PreBalanceProcesses(A3d, ATop, ABottom);
  SplitCElements;
  SplitElements;
  SplitTracers(TracerBased);
  SplitProcesses;
  ApplyRepaints;
  PostBalanceProcesses(A3d, ATop, ABottom, AutoBurialFluxes);
  GenerateIndexes;
end;

end.
