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;
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;
      tracers[length(tracers)-1].isPositive:=0;

      //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;
      tracers[length(tracers)-1].isPositive:=0;

      //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';

{      //adding top tracer
      setLength(tracers,length(tracers)+1);
      initErgTracer(tracers[length(tracers)-1]);
      tracers[length(tracers)-1].name:='aged_'+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].vertLoc:=2;
      tracers[length(tracers)-1].outputUnit:='d*mmol/m**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;

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

procedure SplitElements;
var i: Integer;
begin
  setLength(newCElements,length(cElements));
  for i:=0 to length(cElements)-1 do
  begin
    newCElements[i]:=cElements[i];
    setLength(elements,length(elements)+1);
    elements[length(elements)-1] := elements[cElements[i].myElementNum];
    elements[length(elements)-1].name:=cElements[i].color+'_'+elements[length(elements)-1].name;
    elements[length(elements)-1].description:=cElements[i].description;
    newCElements[i].element:=elements[length(elements)-1].name;
    newCElements[i].myElementNum:=length(elements)-1;
  end;
end;

procedure SplitTracers(tracerBased: Boolean);
var i: Integer;
    t, t_old_max: Integer;
    c: Integer;
    r, rt, rtr: Integer;
    myk,k: Integer;
    gefunden: Boolean;
    isContained: Boolean;
begin
  DecimalSeparator:='.';
  t_old_max:=length(tracers)-1;
  for i:=0 to length(cElements)-1 do
  begin
    for t:=0 to t_old_max do
    begin
      if tracers[t].isActive>0 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].vertLoc:=tracers[t].vertLoc;
          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;

          if tracerBased then
          begin
            tracers[length(tracers)-1].name:=tracers[t].name+'_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';

            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].description:=cElements[i].description+' contained in '+tracers[t].description;

            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;
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.
// Repainting processes are considered as burial of the old colored element
// and reset the age to zero.
var
  c,p, i,o,t,e, r: Integer;
  bal3dI, balTopI, balBottomI: Real;
  bal3dO, balTopO, balBottomO: Real;
begin
  SetLength(A3d     ,length(cElements),length(processes));
  SetLength(ATop    ,length(cElements),length(processes));
  SetLength(ABottom ,length(cElements),length(processes));

  for c:=0 to length(cElements)-1 do
    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 tracers[t].isActive>0 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 tracers[t].isActive>0 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;

procedure SplitProcesses;
// Each process that has colored tracer output will be duplicated.
// It's rate 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 flat tracers this is bottomcellheight/topcellheight to get the units correct
//       [...] means max(0.0,min(1.0,...))
// input needs to be done in another way
//
// However, processes that have colored input only will be duplicated as well. (Burial of element)
// They have no input or output but can contribute to fluxes

const epsilon='0.00000000001';
var ce: Integer;
    p, pp, pmax: Integer;
    i, o: Integer;
    c: Integer;
    t, tt: Integer;
    contained: Boolean;
    ElInputs, cElInputs: TStringList;
    ContentString: String;
begin
  DecimalSeparator:='.';
  ElInputs:=TStringList.create;
  cElInputs:=TStringList.create;
  pmax:=length(processes)-1;
  for ce:=0 to length(cElements)-1 do //loop over colored elements
  begin
    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
          if tracers[processes[p].output[o].myTracerNum].contents[c].myElementNum = cElements[ce].myElementNum then
            if tracers[processes[p].output[o].myTracerNum].isActive>0 then
              contained:=true;
      //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
          if tracers[processes[p].input[o].myTracerNum].contents[c].myElementNum = cElements[ce].myElementNum then
            if tracers[processes[p].input[o].myTracerNum].isActive>0 then
              contained:=true;
      //if contained, the process will be duplicated;
      if contained then
      begin
        setLength(processes,length(processes)+1);
        pp:=length(processes)-1;
        InitErgProcess(processes[pp]);
        processes[pp]:=processes[p];
        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].input,0); //only in output, colored tracers have fixed ratio
        setLength(processes[pp].repaint,0);
        setLength(processes[pp].output,length(processes[p].output));
        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
          if tracers[processes[p].input[i].myTracerNum].isActive>0 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;
        if cElInputs.Count=0 then // no input of colored elements -> deleting process
          setLength(processes,length(processes)-1)
        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;
      end;
    end;
  end;
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
    processes[p]:=processes[p+1];
  setLength(processes,length(processes)-1);
end;

procedure ApplyRepaint(pp, r, e: Integer);
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
    //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
            processes[p].name:='';
    //delete marked subprocesses
    for p:=length(processes)-1 downto 0 do
      if processes[p].name='' then
        DeleteProcess(p);
    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 nessecary, add colored elements to the standard process 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
  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
      //add colored elements to the standard process output, i.e., paint all of the element which passes the process
      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;
      //now 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;
      //nothing happens to the age, that's fine.
    end;
  end
  else
  begin //oldColor <> 'all'/'none'
    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
            if lowercase(processes[pp].repaint[r].newColor)='none' then
            begin
              processes[p].name:='';  //mark for deleting
            end
            else
            begin
              for o:=0 to length(processes[p].output)-1 do
                if pos(lowercase(colorstring),lowercase(processes[p].output[o].tracer))>0 then
                begin
                  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 does not exist -> deleting process
                    processes[p].name:='';
                end;
            end;
        end;
    for p:=length(processes)-1 downto 0 do
      if processes[p].name='' then
        DeleteProcess(p);
  end;
end;

procedure ApplyRepaints;
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);
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;
begin
  DecimalSeparator:='.';
  for c:=0 to length(cElements)-1 do
  begin
    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;
    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;
                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+')';
                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;
                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+')';
                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;
                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+')';
                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;
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.
