unit erg_vector;

interface

procedure SplitVectorTracers(generatePropagationProcesses: Boolean; advancedPropagation: Boolean; ageClassPropagation: Boolean=false);

implementation

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

type TDimensionRec=record
  name: String;
  value: Integer;
end;

var DimArray: Array of TDimensionRec;

procedure CreateDimArray;
var
  t: Integer;
begin
  SetLength(DimArray,0);
  for t:=0 to length(tracers)-1 do
    if tracers[t].dimension>0 then
    begin
      setLength(DimArray, length(DimArray)+1);
      DimArray[length(DimArray)-1].name:='$'+tracers[t].name;
      DimArray[length(DimArray)-1].value:=tracers[t].dimension;
    end;
end;

procedure MakeExplicit(name: String; var s: String);
var
  idx: Array of Integer;
  tempname, newS, s1: String;
  i, d: Integer;
  oldTempname: String;
begin
  tempname:=name;
  //get the DimArray indexes from the name in given order
  setLength(idx,0);
  oldTempname:=Tempname;
  while pos('$',tempname)>0 do
  begin
    for i:=0 to length(DimArray)-1 do
    begin
      if lowercase(copy(tempname,pos('$',tempname),length(DimArray[i].name)))
           =lowercase(DimArray[i].name) then
      begin
        tempname:=copy(tempname,1,pos('$',tempname)-1)+'1'+
                     copy(tempname,pos('$',tempname)+length(DimArray[i].name),length(tempname));
        setLength(idx,length(idx)+1);
        idx[length(idx)-1]:=i;
      end;
    end;
    if tempname=oldTempname then
    begin
      ShowMessage('Error: '+tempname+' could not be split into vector components.');
      tempname:='';
    end;
    oldTempname:=tempname;
  end;
  //now split all appearances of "$something" at the right hand side that also appeared on the left
  //last index changes first
  for i:=0 to length(idx)-1 do
  begin
    newS:='';
    if pos(lowercase(DimArray[idx[i]].name),lowercase(s))>0 then
    begin
      s1:=SemiItem(s);
      while s1<>'' do
      begin
        for d:=1 to DimArray[idx[i]].value do
        begin
          newS:=newS+StringReplace(s1,DimArray[idx[i]].name,IntToStr(d),[rfReplaceAll, rfIgnoreCase])+';';
        end;
        s1:=SemiItem(s);
      end;
    end
    else   //index appears at the left only
    begin
      newS:=s;
    end;
    if copy(newS,length(newS),1)=';' then
      s:=copy(newS,1,length(newS)-1) //remove final ";"
    else
      s:=newS;
  end;
  //now get the dimensions that remain at the right hand side
  setLength(idx,0);
  for i:=0 to length(DimArray)-1 do
  begin
    if pos(DimArray[i].name,s)>0 then
    begin
      setLength(idx,length(idx)+1);
      idx[length(idx)-1]:=i;
    end;
  end;
  //now perform summation over all these indexes
  newS:='';
  for i:=0 to length(idx)-1 do
  begin
    s1:=SemiItem(s);
    while s1<>'' do
    begin
      if pos(lowercase(DimArray[idx[i]].name),lowercase(s1))>0 then
      begin
        for d:=1 to DimArray[idx[i]].value do
        begin
          newS:=newS+'('+StringReplace(s1,DimArray[idx[i]].name,IntToStr(d),[rfReplaceAll, rfIgnoreCase])+')+';
        end;
        newS:=copy(newS,1,length(newS)-1)+';';
      end
      else
        newS:=newS+IntToStr(DimArray[idx[i]].value)+'*('+s1+');';
      s1:=SemiItem(s);
    end;
    if copy(newS,length(newS),1)=';' then
      s:=copy(newS,1,length(newS)-1) //remove final ";"
    else
      s:=newS;
  end;
end;

procedure MakeConstantsExplicit;
var
  c: Integer;
begin
  for c:=0 to length(constants)-1 do
  begin
    //For vector constants, we will check if all vector elements are the same.
    //This is the case if no ";" splits different entries.
    //If all are the same, we will set myIsMissingInDocumentation to true.
    //This will later be changed in the first element only,
    //so only this one appears in the documentation in which we filter by
    //<auxiliaries isMissingInDocumentation=0>
    //</auxiliaries>
    constants[c].myIsMissingInDocumentation:=false;
    if pos('$',constants[c].name)>0 then constants[c].myIsMissingInDocumentation:=true;
    if pos(';',constants[c].valueString)>0 then constants[c].myIsMissingInDocumentation:=false;

    //now split the right-hand sides into vector elements, separated by ";"
    //e.g.:
    //  a_$t_det = b_$t_det + 1
    //will become
    //  a_$t_det = b_1 + 1; b_2 + 1; ...; b_6 + 1
    MakeExplicit(constants[c].name,constants[c].valueString);
  end;
end;

procedure MakeAuxiliariesExplicit;
var
  a, i: Integer;
begin
  for a:=0 to length(auxiliaries)-1 do
  begin
    //For vector auxiliaries, we will check if all vector elements are the same.
    //This is the case if no ";" splits different entries.
    //If all are the same, we will set myIsMissingInDocumentation to true.
    //This will later be changed in the first element only,
    //so only this one appears in the documentation in which we filter by
    //<auxiliaries isMissingInDocumentation=0>
    //</auxiliaries>
    auxiliaries[a].myIsMissingInDocumentation:=false;
    if pos('$',auxiliaries[a].name)>0 then auxiliaries[a].myIsMissingInDocumentation:=true;
    if pos(';',auxiliaries[a].formula)>0 then auxiliaries[a].myIsMissingInDocumentation:=false;
    for i:=1 to 9 do
      if pos(';',auxiliaries[a].temp[i])>0 then auxiliaries[a].myIsMissingInDocumentation:=false;

    //now split the right-hand sides into vector elements, separated by ";"
    //e.g.:
    //  a_$t_det = b_$t_det + 1
    //will become
    //  a_$t_det = b_1 + 1; b_2 + 1; ...; b_6 + 1
    for i:=1 to 9 do
      MakeExplicit(auxiliaries[a].name,auxiliaries[a].temp[i]);
    MakeExplicit(auxiliaries[a].name,auxiliaries[a].formula);
  end;
end;

procedure MakeProcessesExplicit;
var
  p: Integer;
begin
  for p:=0 to length(processes)-1 do
  begin
    //For vector processes, we will check if all vector elements are the same.
    //This is the case if no ";" splits different entries.
    //If all are the same, we will set myIsMissingInDocumentation to true.
    //This will later be changed in the first element only,
    //so only this one appears in the documentation in which we filter by
    //<auxiliaries isMissingInDocumentation=0>
    //</auxiliaries>
    processes[p].myIsMissingInDocumentation:=false;
    if pos('$',processes[p].name)>0 then processes[p].myIsMissingInDocumentation:=true;
    if pos(';',processes[p].turnover)>0 then processes[p].myIsMissingInDocumentation:=false;

    //now split the right-hand sides into vector elements, separated by ";"
    //e.g.:
    //  a_$t_det = b_$t_det + 1
    //will become
    //  a_$t_det = b_1 + 1; b_2 + 1; ...; b_6 + 1
    MakeExplicit(processes[p].name,processes[p].turnover);
  end;
end;

procedure SplitConstants;
var
  allAreSplit: Boolean;
  firstToSplit: Integer;
  i, j: Integer;
  idx: Array of Integer;
  currentIdx: Array of Integer;
  tempname: String;
  product: Integer;
  s: String;
begin
  allAreSplit:=true;
  for i:=0 to length(constants)-1 do
    if pos('$',constants[i].name)>0 then
      if allAreSplit then
      begin
        allAreSplit:=false;
        firstToSplit:=i;
      end;
  while allAreSplit=false do
  begin
    // generate list of dimensions contained
    setLength(idx,0);
    tempname:=constants[firstToSplit].name;
    while pos('$',tempname)>0 do
    begin
      for i:=0 to length(DimArray)-1 do
      begin
        if lowercase(copy(tempname,pos('$',tempname),length(DimArray[i].name)))
             =lowercase(DimArray[i].name) then
        begin
          tempname:=copy(tempname,1,pos('$',tempname)-1)+'1'+
                       copy(tempname,pos('$',tempname)+length(DimArray[i].name),length(tempname));
          setLength(idx,length(idx)+1);
          idx[length(idx)-1]:=i;
        end;
      end;
    end;

    //now calculate product of all dimensions
    product:=1;
    for i:=0 to length(idx)-1 do
      product:=product*DimArray[idx[i]].value;

    //multiply the existing constant and make some space for its siblings
    //by moving the other constants down in the list
    setLength(constants,length(constants)+product-1);
    for i:=length(constants)-1 downto firstToSplit+product do
      copyErgConstant(constants[i-product+1],constants[i]);
    for i:=firstToSplit+1 to firstToSplit+product-1 do
      copyErgConstant(constants[firstToSplit],constants[i]);

    //if isMissingInDocumentation was set, which means later siblings are hidden
    //from the documentation, do two things:
    //- change it back in the first one, so this one appears in the documentation
    //- change that one's description so it states there are these other siblings
    if constants[firstToSplit].myIsMissingInDocumentation=true then
    begin
      constants[firstToSplit].myDescriptionInDocumentation:=constants[firstToSplit].description;
      for i:=0 to length(idx)-1 do
        constants[firstToSplit].myDescriptionInDocumentation:=
          constants[firstToSplit].myDescriptionInDocumentation + ' (index 1 to '+IntToStr(DimArray[idx[i]].value)+')';
      constants[firstToSplit].myIsMissingInDocumentation:=false;
    end;

    //split the ";"-separated values
    s:=constants[firstToSplit].valueString;
    for i:=firstToSplit to firstToSplit+product-1 do
      constants[i].valueString:=SemiItem(s);

    //replace "$"-expressions in name by indexes, last index changes first
    setLength(currentIdx,length(idx));
    for i:=0 to length(currentIdx)-1 do
      currentIdx[i]:=1;
    for i:=firstToSplit to firstToSplit+product-1 do
    begin
      for j:=0 to length(idx)-1 do
        constants[i].name:=StringReplace(constants[i].name,DimArray[idx[j]].name,IntToStr(currentIdx[j]),[rfReplaceAll, rfIgnoreCase]);
      currentIdx[length(currentIdx)-1]:=currentIdx[length(currentIdx)-1]+1;
      for j:=length(currentIdx)-1 downto 1 do
        if currentIdx[j]>DimArray[idx[j]].value then
        begin
          currentIdx[j-1]:=currentIdx[j-1]+1;
          currentIdx[j]:=1;
        end;
    end;

    allAreSplit:=true;
    for i:=0 to length(constants)-1 do
      if pos('$',constants[i].name)>0 then
        if allAreSplit then
        begin
          allAreSplit:=false;
          firstToSplit:=i;
        end;
  end;
end;

procedure SplitAuxiliaries;
var
  allAreSplit: Boolean;
  firstToSplit: Integer;
  i, j, k: Integer;
  idx: Array of Integer;
  currentIdx: Array of Integer;
  tempname: String;
  product: Integer;
  s: String;
begin
  allAreSplit:=true;
  for i:=0 to length(Auxiliaries)-1 do
    if pos('$',Auxiliaries[i].name)>0 then
      if allAreSplit then
      begin
        allAreSplit:=false;
        firstToSplit:=i;
      end;
  while allAreSplit=false do
  begin
    // generate list of dimensions contained
    setLength(idx,0);
    tempname:=Auxiliaries[firstToSplit].name;
    while pos('$',tempname)>0 do
    begin
      for i:=0 to length(DimArray)-1 do
      begin
        if lowercase(copy(tempname,pos('$',tempname),length(DimArray[i].name)))
             =lowercase(DimArray[i].name) then
        begin
          tempname:=copy(tempname,1,pos('$',tempname)-1)+'1'+
                       copy(tempname,pos('$',tempname)+length(DimArray[i].name),length(tempname));
          setLength(idx,length(idx)+1);
          idx[length(idx)-1]:=i;
        end;
      end;
    end;

    //now calculate product of all dimensions
    product:=1;
    for i:=0 to length(idx)-1 do
      product:=product*DimArray[idx[i]].value;

    //multiply the existing Auxiliary and make some space for its siblings
    //by moving the other auxiliaries down in the list
    setLength(Auxiliaries,length(Auxiliaries)+product-1);
    for i:=length(Auxiliaries)-1 downto firstToSplit+product do
      copyErgAuxiliary(Auxiliaries[i-product+1],Auxiliaries[i]);
    for i:=firstToSplit+1 to firstToSplit+product-1 do
      copyErgAuxiliary(Auxiliaries[firstToSplit],Auxiliaries[i]);

    //if isMissingInDocumentation was set, which means later siblings are hidden
    //from the documentation, do two things:
    //- change it back in the first one, so this one appears in the documentation
    //- change that one's description so it states there are these other siblings
    if auxiliaries[firstToSplit].myIsMissingInDocumentation=true then
    begin
      auxiliaries[firstToSplit].myDescriptionInDocumentation:=auxiliaries[firstToSplit].description;
      for i:=0 to length(idx)-1 do
        auxiliaries[firstToSplit].myDescriptionInDocumentation:=
          auxiliaries[firstToSplit].myDescriptionInDocumentation + ' (index 1 to '+IntToStr(DimArray[idx[i]].value)+')';
      auxiliaries[firstToSplit].myIsMissingInDocumentation:=false;
    end;

    //split the ";"-separated values
    s:=Auxiliaries[firstToSplit].formula;
    for i:=firstToSplit to firstToSplit+product-1 do
      Auxiliaries[i].formula:=SemiItem(s);

    for k:=1 to 9 do
    begin
      s:=Auxiliaries[firstToSplit].temp[k];
      for i:=firstToSplit to firstToSplit+product-1 do
        Auxiliaries[i].temp[k]:=SemiItem(s);
    end;

    //replace "$"-expressions in name by indexes, last index changes first
    setLength(currentIdx,length(idx));
    for i:=0 to length(currentIdx)-1 do
      currentIdx[i]:=1;
    for i:=firstToSplit to firstToSplit+product-1 do
    begin
      for j:=0 to length(idx)-1 do
        Auxiliaries[i].name:=StringReplace(Auxiliaries[i].name,DimArray[idx[j]].name,IntToStr(currentIdx[j]),[rfReplaceAll, rfIgnoreCase]);
      currentIdx[length(currentIdx)-1]:=currentIdx[length(currentIdx)-1]+1;
      for j:=length(currentIdx)-1 downto 1 do
        if currentIdx[j]>DimArray[idx[j]].value then
        begin
          currentIdx[j-1]:=currentIdx[j-1]+1;
          currentIdx[j]:=1;
        end;
    end;

    allAreSplit:=true;
    for i:=0 to length(Auxiliaries)-1 do
      if pos('$',Auxiliaries[i].name)>0 then
        if allAreSplit then
        begin
          allAreSplit:=false;
          firstToSplit:=i;
        end;
  end;
end;

procedure SplitProcesses;
var
  allAreSplit: Boolean;
  firstToSplit: Integer;
  i, j, k: Integer;
  idx: Array of Integer;
  currentIdx: Array of Integer;
  tempname: String;
  product: Integer;
  s: String;
begin
  allAreSplit:=true;
  for i:=0 to length(Processes)-1 do
    if pos('$',Processes[i].name)>0 then
      if allAreSplit then
      begin
        allAreSplit:=false;
        firstToSplit:=i;
      end;
  while allAreSplit=false do
  begin
    // generate list of dimensions contained
    setLength(idx,0);
    tempname:=Processes[firstToSplit].name;
    while pos('$',tempname)>0 do
    begin
      for i:=0 to length(DimArray)-1 do
      begin
        if lowercase(copy(tempname,pos('$',tempname),length(DimArray[i].name)))
             =lowercase(DimArray[i].name) then
        begin
          tempname:=copy(tempname,1,pos('$',tempname)-1)+'1'+
                       copy(tempname,pos('$',tempname)+length(DimArray[i].name),length(tempname));
          setLength(idx,length(idx)+1);
          idx[length(idx)-1]:=i;
        end;
      end;
    end;

    //now calculate product of all dimensions
    product:=1;
    for i:=0 to length(idx)-1 do
      product:=product*DimArray[idx[i]].value;

    //multiply the existing Process and make some space for its siblings
    //by moving the other processes down in the list
    setLength(Processes,length(Processes)+product-1);
    for i:=length(Processes)-1 downto firstToSplit+product do
      copyErgProcess(Processes[i-product+1],Processes[i]);
    for i:=firstToSplit+1 to firstToSplit+product-1 do
      copyErgProcess(Processes[firstToSplit],Processes[i]);

    //if isMissingInDocumentation was set, which means later siblings are hidden
    //from the documentation, do two things:
    //- change it back in the first one, so this one appears in the documentation
    //- change that one's description so it states there are these other siblings
    if processes[firstToSplit].myIsMissingInDocumentation=true then
    begin
      processes[firstToSplit].myDescriptionInDocumentation:=processes[firstToSplit].description;
      for i:=0 to length(idx)-1 do
        processes[firstToSplit].myDescriptionInDocumentation:=
          processes[firstToSplit].myDescriptionInDocumentation + ' (index 1 to '+IntToStr(DimArray[idx[i]].value)+')';
      processes[firstToSplit].myIsMissingInDocumentation:=false;
    end;

    //split the ";"-separated values
    s:=Processes[firstToSplit].turnover;
    for i:=firstToSplit to firstToSplit+product-1 do
      Processes[i].turnover:=SemiItem(s);

    //replace "$"-expressions in name and input/output by indexes, last index changes first
    setLength(currentIdx,length(idx));
    for i:=0 to length(currentIdx)-1 do
      currentIdx[i]:=1;
    for i:=firstToSplit to firstToSplit+product-1 do
    begin
      for j:=0 to length(idx)-1 do
      begin
        Processes[i].name:=StringReplace(Processes[i].name,DimArray[idx[j]].name,IntToStr(currentIdx[j]),[rfReplaceAll, rfIgnoreCase]);
        for k:=0 to length(processes[i].input)-1 do
        begin
          Processes[i].input[k].tracer:=StringReplace(Processes[i].input[k].tracer,DimArray[idx[j]].name,IntToStr(currentIdx[j]),[rfReplaceAll, rfIgnoreCase]);
          Processes[i].input[k].amount:=StringReplace(Processes[i].input[k].amount,DimArray[idx[j]].name,IntToStr(currentIdx[j]),[rfReplaceAll, rfIgnoreCase]);
        end;
        for k:=0 to length(processes[i].output)-1 do
        begin
          Processes[i].output[k].tracer:=StringReplace(Processes[i].output[k].tracer,DimArray[idx[j]].name,IntToStr(currentIdx[j]),[rfReplaceAll, rfIgnoreCase]);
          Processes[i].output[k].amount:=StringReplace(Processes[i].output[k].amount,DimArray[idx[j]].name,IntToStr(currentIdx[j]),[rfReplaceAll, rfIgnoreCase]);
        end;
        for k:=0 to length(processes[i].limitations)-1 do
        begin
          //needs to be done to apply posdef scheme
        end;
      end;
      currentIdx[length(currentIdx)-1]:=currentIdx[length(currentIdx)-1]+1;
      for j:=length(currentIdx)-1 downto 1 do
        if currentIdx[j]>DimArray[idx[j]].value then
        begin
          currentIdx[j-1]:=currentIdx[j-1]+1;
          currentIdx[j]:=1;
        end;
    end;

    allAreSplit:=true;
    for i:=0 to length(Processes)-1 do
      if pos('$',Processes[i].name)>0 then
        if allAreSplit then
        begin
          allAreSplit:=false;
          firstToSplit:=i;
        end;
  end;
end;

procedure SplitTracers(generatePropagationProcesses: Boolean; advancedPropagation: Boolean; ageClassPropagation: Boolean=false);
var
  allAreSplit, found: Boolean;
  firstToSplit: Integer;
  i, j, k: Integer;
  idx: Array of Integer;
  currentIdx: Array of Integer;
  tempname: String;
  product: Integer;
  s: String;
  oldTracerName, MassLimitString, lowerMassLimitString, upperMassLimitString: String;
  g, h, meanMassStationary, cfraction, bs, meanMass, b, c2: String;
begin
  decimalSeparator:='.';
  allAreSplit:=true;
  for i:=0 to length(Tracers)-1 do
    if tracers[i].dimension>0 then
      if allAreSplit then
      begin
        allAreSplit:=false;
        firstToSplit:=i;
      end;
  while allAreSplit=false do
  begin

    //now calculate product of all dimensions
    product:=tracers[FirstToSplit].dimension;
    if tracers[FirstToSplit].massLimits <> '' then
      product:=product*2;

    //multiply the existing Tracer
    setLength(Tracers,length(Tracers)+product-1);
    for i:=length(Tracers)-1 downto firstToSplit+product do
      copyErgTracer(Tracers[i-product+1],Tracers[i]);
    for i:=firstToSplit+1 to firstToSplit+product-1 do
      copyErgTracer(Tracers[firstToSplit],Tracers[i]);

    //replace "$"-expressions in name by indexes, last index changes first
    for i:=firstToSplit to firstToSplit+product-1 do
    begin
      j:=i-firstToSplit;
      if tracers[i].massLimits <> '' then j:=j div 2;
      j:=j+1;
      Tracers[i].name:=Tracers[i].name+'_'+IntToStr(j);
      Tracers[i].description:=Tracers[i].description+'_'+IntToStr(j);
      tracers[i].dimension:=0;
    end;
    //split initValue: mass1;abd1;mass2;abd2;...
    s:=tracers[firstToSplit].initValue;
    for i:=firstToSplit to firstToSplit+product-1 do
      tracers[i].initValue:=SemiItem(s);
    if tracers[FirstToSplit].massLimits<>'' then  //stage-resolved vector tracer
    begin
      for i:=firstToSplit to firstToSplit+product-1 do    //create mass and abundance tracers
      begin
        j:=i-firstToSplit;
        j:=j mod 2;
        if j=0 then
        begin
          Tracers[i].name:=Tracers[i].name+'_mass';
          Tracers[i].description:=Tracers[i].description+'_mass';
        end
        else
        begin
          Tracers[i].name:=Tracers[i].name+'_abd';
          Tracers[i].description:=Tracers[i].description+'_abd';
          setLength(tracers[i].contents,0);
        end;
      end;
  {    for i:=1 to (product div 2) do                    //define average mass auxiliary
      begin
        oldTracerName:=copy(tracers[firstToSplit].name,1,length(tracers[firstToSplit].name)-7);
        SetLength(auxiliaries,length(auxiliaries)+1);
        InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
        auxiliaries[length(auxiliaries)-1].name:='avg_mass_'+OldTracerName+'_'+IntToStr(i);
        auxiliaries[length(auxiliaries)-1].formula:='theta('+OldTracerName+'_'+IntToStr(i)+'_abd)*'+OldTracerName+'_'+IntToStr(i)+'_mass / max(1.0e-34,'+OldTracerName+'_'+IntToStr(i)+'_abd)';
        auxiliaries[length(auxiliaries)-1].vertLoc:=tracers[firstToSplit].vertLoc;
        auxiliaries[length(auxiliaries)-1].description:='average mass of '+tracers[firstToSplit].description+' class '+IntToStr(i)+' [mmol]';
      end;     }
      MassLimitString:=tracers[firstToSplit].massLimits;
      lowerMassLimitString:=SemiItem(MassLimitString);
      if generatePropagationProcesses then
      begin
        //first, generate an auxiliary variable for the total time step to avoid conflicts with positive-definite scheme which uses partial steps
        found:=false;
        for i:=0 to length(auxiliaries)-1 do
          if auxiliaries[i].name='cgt_timestep_total' then
            found:=true;
        if not found then
        begin
          setLength(auxiliaries,length(auxiliaries)+1);
          initErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
          auxiliaries[length(auxiliaries)-1].name:='cgt_timestep_total';
          auxiliaries[length(auxiliaries)-1].formula:='cgt_timestep';
        end;
        for i:=1 to (product div 2) do
        begin
         upperMassLimitString:=SemiItem(MassLimitString);
         oldTracerName:=copy(tracers[firstToSplit].name,1,length(tracers[firstToSplit].name)-7);
         //generate constants: tracer_$tracer_min_mass, tracer_$tracer_max_mass
         SetLength(constants,length(constants)+1);
         InitErgConstant(constants[length(constants)-1]);
         constants[length(constants)-1].name:=oldTracerName+'_'+IntToStr(i)+'_min_mass';
         constants[length(constants)-1].valueString:=LowerMassLimitString;
         constants[length(constants)-1].description:='lower limit for individual mass of '+oldTracerName+' of class '+IntToStr(i);
         if upperMassLimitString<>'' then
         begin
           SetLength(constants,length(constants)+1);
           InitErgConstant(constants[length(constants)-1]);
           constants[length(constants)-1].name:=oldTracerName+'_'+IntToStr(i)+'_max_mass';
           constants[length(constants)-1].valueString:=UpperMassLimitString;
           constants[length(constants)-1].description:='upper limit for individual mass of '+oldTracerName+' of class '+IntToStr(i);
         end;
         if i <=(product div 2)-1 then //stage propagation processes
         begin
          if advancedPropagation then
          begin
            g:='max(1.0e-10,propagation_'+oldTracerName+'_growth_'+IntToStr(i)+')';
            h:='max(0.0,propagation_'+oldTracerName+'_mortality_'+IntToStr(i)+')';
            meanMassStationary:='theta('+h+'-1.0e-5)*theta(power('+g+'-'+h+',2.0)-1.0e-10)*'+
                                     'max('+lowerMassLimitString+',min('+upperMassLimitString+','+h+'/('+g+'-'+h+')*(power('+upperMassLimitString+',1.0-'+h+'/'+g+')-power('+lowerMassLimitString+',1.0-'+h+'/'+g+'))/'
                                    +'(power('+lowerMassLimitString+',max(1.0e-5,0.0-'+h+'/'+g+'))-power('+upperMassLimitString+',max(1.0e-5,0.0-'+h+'/'+g+')))))'
                                 +'+(1.0-theta('+h+'-1.0e-5))*theta(power('+g+'-'+h+',2.0)-1.0e-10)*'+
                                    '('+upperMassLimitString+'-'+lowerMassLimitString+')/log('+upperMassLimitString+'/('+lowerMassLimitString+'))'
                                 +'+(1.0-theta(power('+g+'-'+h+',2.0)-1.0e-10))*'+
                                    'log('+upperMassLimitString+'/('+lowerMassLimitString+'))/('+upperMassLimitString+'-'+lowerMassLimitString+')';

            SetLength(processes,length(processes)+1);
            InitErgProcess(processes[length(processes)-1]);
            processes[length(processes)-1].name:='propagation_'+oldTracerName+'_meanmassstationary_'+IntToStr(i)+'_to_'+IntToStr(i+1);
            processes[length(processes)-1].description:='relative growth rate of '+tracers[firstToSplit].description+' in mass class '+IntToStr(i)+' [1/day]';
            processes[length(processes)-1].vertLoc:=tracers[firstToSplit].vertLoc;
            setLength(processes[length(processes)-1].input,0);
            setLength(processes[length(processes)-1].output,0);
            processes[length(processes)-1].turnover:=MeanMassStationary;
            processes[length(processes)-1].processType:='propagation';
            meanMassStationary:=processes[length(processes)-1].name;

            cfraction:='theta('+h+'-1.0e-5)*'+
                           'max(0.0,power('+upperMassLimitString+',0.0-1.0-'+h+'/'+g+')*'+h+'/'+g+'*('+upperMassLimitString+'-'+lowerMassLimitString+')/'
                           +'(power('+lowerMassLimitString+',0.0-'+h+'/'+g+')-power('+upperMassLimitString+',0.0-'+h+'/'+g+')))'
                       +'+(1-theta('+h+'-1.0e-5))*'+
                           '(1.0-'+lowerMassLimitString+'/('+upperMassLimitString+'))/log('+upperMassLimitString+'/('+lowerMassLimitString+'))';
            BS:='max(0.0,'+g+'*('+meanMassStationary+')-'+g+'*('+cfraction+')*'+upperMassLimitString+'*('+upperMassLimitString+'-('+meanMassStationary+'))/'
                +'('+upperMassLimitString+'-'+lowerMassLimitString+'))';
            SetLength(processes,length(processes)+1);
            InitErgProcess(processes[length(processes)-1]);
            processes[length(processes)-1].name:='propagation_'+oldTracerName+'_bs_'+IntToStr(i)+'_to_'+IntToStr(i+1);
            processes[length(processes)-1].description:='relative growth rate of '+tracers[firstToSplit].description+' in mass class '+IntToStr(i)+' [1/day]';
            processes[length(processes)-1].vertLoc:=tracers[firstToSplit].vertLoc;
            setLength(processes[length(processes)-1].input,0);
            setLength(processes[length(processes)-1].output,0);
            processes[length(processes)-1].turnover:=BS;
            processes[length(processes)-1].processType:='propagation';
            BS:=processes[length(processes)-1].name;

            meanMass:='max('+lowerMassLimitString+',min('+upperMassLimitString+','+'theta('+OldTracerName+'_'+IntToStr(i)+'_abd)*'+OldTracerName+'_'+IntToStr(i)+'_mass / max(1.0e-30,'+OldTracerName+'_'+IntToStr(i)+'_abd)))';
            B:='theta('+meanMassStationary+'-'+meanMass+')*('+g+'*'+lowerMassLimitString+'+(('+BS+')-'+g+'*'+lowerMassLimitString+')*('+meanMass+'-'+lowerMassLimitString+')/max(1.0e-10,('+meanMassStationary+')-'+lowerMassLimitString+'))'
               +'+theta('+meanMass+'-'+MeanMassStationary+')*(('+BS+')*(1.0-('+meanMass+'-'+meanMassStationary+')/max(1.0e-10,'+upperMassLimitString+'-'+meanMassStationary+')))';
            SetLength(processes,length(processes)+1);
            InitErgProcess(processes[length(processes)-1]);
            processes[length(processes)-1].name:='propagation_'+oldTracerName+'_b_'+IntToStr(i)+'_to_'+IntToStr(i+1);
            processes[length(processes)-1].description:='relative growth rate of '+tracers[firstToSplit].description+' in mass class '+IntToStr(i)+' [1/day]';
            processes[length(processes)-1].vertLoc:=tracers[firstToSplit].vertLoc;
            setLength(processes[length(processes)-1].input,0);
            setLength(processes[length(processes)-1].output,0);
            processes[length(processes)-1].turnover:=B;
            processes[length(processes)-1].processType:='propagation';
            B:=processes[length(processes)-1].name;

          end;
          SetLength(processes,length(processes)+1);
          InitErgProcess(processes[length(processes)-1]);
          processes[length(processes)-1].name:='propagation_'+oldTracerName+'_mass_'+IntToStr(i)+'_to_'+IntToStr(i+1);
          processes[length(processes)-1].description:='Propagation of '+tracers[firstToSplit].description+' mass from mass class '+IntToStr(i)+' to '+IntToStr(i+1);
          processes[length(processes)-1].vertLoc:=tracers[firstToSplit].vertLoc;
          setLength(processes[length(processes)-1].input,1);
          processes[length(processes)-1].input[0].tracer:=oldTracerName+'_'+IntToStr(i)+'_mass';
          processes[length(processes)-1].input[0].amount:='1.0';
          setLength(processes[length(processes)-1].output,1);
          processes[length(processes)-1].output[0].tracer:=oldTracerName+'_'+IntToStr(i+1)+'_mass';
          processes[length(processes)-1].output[0].amount:='1.0';
          processes[length(processes)-1].processType:='propagation';
          if ageClassPropagation=true then
          begin
            processes[length(processes)-1].turnover:=oldTracerName+'_'+IntToStr(i)+'_mass'+
               '*(theta(cgt_timestep*1.5-(cgt_dayofyear-1+cgt_hour/24.0))-theta(cgt_timestep*0.5-(cgt_dayofyear-1+cgt_hour/24.0)))/cgt_timestep_total';
          end
          else if advancedPropagation=false then
          begin
            processes[length(processes)-1].turnover:=oldTracerName+'_'+IntToStr(i)+'_mass*theta('+
                                               'theta('+OldTracerName+'_'+IntToStr(i)+'_abd)*'+OldTracerName+'_'+IntToStr(i)+'_mass / max(1.0e-30,'+OldTracerName+'_'+IntToStr(i)+'_abd)'
                                               +'-('+upperMassLimitString+'))/cgt_timestep_total';
          end
          else
          begin
            //c2=OldTracerName+'_'+IntToStr(i)+'_abd*(g*meanMass-B)/(g*m2*(m2-meanMass))*theta(meanMass-m1)
            c2:=OldTracerName+'_'+IntToStr(i)+'_abd*('+g+'*'+meanMass+'-('+B+'))/'+
               'max(1.0e-10,'+g+'*'+upperMassLimitString+'*('+upperMassLimitString+'-('+meanMass+')))*theta('+meanMass+'-'+lowerMassLimitString+')';
            processes[length(processes)-1].turnover:='(1-theta('+
                                               'theta('+OldTracerName+'_'+IntToStr(i)+'_abd)*'+OldTracerName+'_'+IntToStr(i)+'_mass / max(1.0e-30,'+OldTracerName+'_'+IntToStr(i)+'_abd)'
                                               +'-('+upperMassLimitString+')))*max(0.0,'+c2+'*'+g+'*'+upperMassLimitString+')'
                                               +'+'+oldTracerName+'_'+IntToStr(i)+'_mass*theta('+
                                              'theta('+OldTracerName+'_'+IntToStr(i)+'_abd)*'+OldTracerName+'_'+IntToStr(i)+'_mass / max(1.0e-30,'+OldTracerName+'_'+IntToStr(i)+'_abd)'
                                               +'-('+upperMassLimitString+'))/cgt_timestep_total';
          end;

          SetLength(processes,length(processes)+1);
          InitErgProcess(processes[length(processes)-1]);
          processes[length(processes)-1].name:='propagation_'+oldTracerName+'_abd_'+IntToStr(i)+'_to_'+IntToStr(i+1);
          processes[length(processes)-1].description:='Propagation of '+tracers[firstToSplit].description+' abundance from mass class '+IntToStr(i)+' to '+IntToStr(i+1);
          processes[length(processes)-1].vertLoc:=tracers[firstToSplit].vertLoc;
          setLength(processes[length(processes)-1].input,1);
          processes[length(processes)-1].input[0].tracer:=oldTracerName+'_'+IntToStr(i)+'_abd';
          processes[length(processes)-1].input[0].amount:='1.0';
          setLength(processes[length(processes)-1].output,1);
          processes[length(processes)-1].output[0].tracer:=oldTracerName+'_'+IntToStr(i+1)+'_abd';
          processes[length(processes)-1].output[0].amount:='1.0';
          processes[length(processes)-1].processType:='propagation';
          if ageClassPropagation=true then
          begin
            processes[length(processes)-1].turnover:=oldTracerName+'_'+IntToStr(i)+'_abd'+
               '*(theta(cgt_timestep*1.5-(cgt_dayofyear-1+cgt_hour/24.0))-theta(cgt_timestep*0.5-(cgt_dayofyear-1+cgt_hour/24.0)))/cgt_timestep_total';
          end
          else if advancedPropagation=false then
            processes[length(processes)-1].turnover:=oldTracerName+'_'+IntToStr(i)+'_abd*theta('+
                                               'theta('+OldTracerName+'_'+IntToStr(i)+'_abd)*'+OldTracerName+'_'+IntToStr(i)+'_mass / max(1.0e-30,'+OldTracerName+'_'+IntToStr(i)+'_abd)'
                                               +'-('+upperMassLimitString+'))/cgt_timestep_total'
          else
          begin
            //c2=OldTracerName+'_'+IntToStr(i)+'_abd*(g*meanMass-B)/(g*m2*(m2-meanMass))
            c2:=OldTracerName+'_'+IntToStr(i)+'_abd*('+g+'*'+meanMass+'-('+B+'))/'+
               'max(1.0e-10,'+g+'*'+upperMassLimitString+'*('+upperMassLimitString+'-('+meanMass+')))*theta('+meanMass+'-'+lowerMassLimitString+')';
            processes[length(processes)-1].turnover:='(1-theta('+
                                               'theta('+OldTracerName+'_'+IntToStr(i)+'_abd)*'+OldTracerName+'_'+IntToStr(i)+'_mass / max(1.0e-30,'+OldTracerName+'_'+IntToStr(i)+'_abd)'
                                               +'-('+upperMassLimitString+')))*max(0.0,'+c2+'*'+g+')'
                                               +'+'
                                               +OldTracerName
                                               +'_'
                                               +IntToStr(i)
                                               +'_abd*theta('
                                               +'theta('+OldTracerName+'_'+IntToStr(i)+'_abd)*'+OldTracerName+'_'+IntToStr(i)+'_mass / max(1.0e-30,'+OldTracerName+'_'+IntToStr(i)+'_abd)'
                                               +'-('+upperMassLimitString+'))/cgt_timestep_total';
          end;
         end;
         lowerMassLimitString:=upperMassLimitString;

        end;
      end;
    end;


    allAreSplit:=true;
    for i:=0 to length(Tracers)-1 do
      if Tracers[i].dimension>0 then
        if allAreSplit then
        begin
          allAreSplit:=false;
          firstToSplit:=i;
        end;
  end;
end;

procedure SplitVectorTracers(generatePropagationProcesses: Boolean; advancedPropagation: Boolean; ageClassPropagation: Boolean=false);
begin
  //first create array which stores the dimension lenghts, e.g. $t_det = 6
  CreateDimArray;

  //Split the right-hand sides
  //e.g.:
  //  a_$t_det = b_$t_det + 1
  //will become
  //  a_$t_det = b_1 + 1; b_2 + 1; ...; b_6 + 1
  MakeConstantsExplicit;
  MakeAuxiliariesExplicit;
  MakeProcessesExplicit;

  //Split the left-hand sides
  //e.g.:
  //  a_$t_det = b_1 + 1; b_2 + 1; ...; b_6 + 1
  //will become
  //  a_1 = b_1 + 1
  //  a_2 = b_2 + 1
  //  ...
  //  a_6 = b_6 + 1
  SplitConstants;
  SplitAuxiliaries;
  SplitProcesses;

  SplitTracers(generatePropagationProcesses, advancedPropagation, ageClassPropagation);

  GenerateIndexes;
end;

end.
