unit erg_base;

interface

uses erg_types, classes, Dialogs;

const
  NullLimitationArray: Array[0..0] of TErgLimitation = (());

var
  constants:   Array of tErgConstant;
  elements:    Array of tErgElement;
  tracers:     Array of tErgTracer;
  auxiliaries: Array of tErgAuxiliary;
  limitations: Array of tErgLimitation;
  processes:   Array of tErgProcess;
  cElements:   Array of tErgCElement;
  modelinfos:  tErgModelInfo;
  keywords:    Array of String;
  keywordDescriptions: Array of String;

var ConstantsList, AuxiliariesList, ProcessesList: TStringList;
  maxIterations: Integer;

function LoadConstants(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
procedure SaveConstants(filename: String);
function LoadSingleConstant(var F: Textfile):Integer; //returns 0 if succeeded
procedure SaveSingleConstant(var F: Textfile; i: Integer);
procedure SaveSingleConstantRef(var F: TextFile; i: Integer; c: TErgConstant; AlwaysShowDefault: Boolean=true);

function LoadElements(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
procedure SaveElements(filename: String);
function LoadSingleElement(var F: Textfile):Integer; //returns 0 if succeeded
procedure SaveSingleElement(var F: Textfile; i: Integer);
procedure SaveSingleElementRef(var F: TextFile; i: Integer; c: TErgElement; AlwaysShowDefault: Boolean=true);

function LoadTracers(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
procedure SaveTracers(filename: String);
function LoadSingleTracer(var F: Textfile):Integer; //returns 0 if succeeded
procedure SaveSingleTracer(var F: Textfile; i: Integer);
procedure SaveSingleTracerRef(var F: TextFile; i: Integer; c: TErgTracer; AlwaysShowDefault: Boolean=true);

function LoadAuxiliaries(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
procedure SaveAuxiliaries(filename: String);
function LoadSingleAuxiliary(var F: Textfile):Integer; //returns 0 if succeeded
procedure SaveSingleAuxiliary(var F: Textfile; i: Integer);
procedure SaveSingleAuxiliaryRef(var F: TextFile; i: Integer; c: TErgAuxiliary; AlwaysShowDefault: Boolean=true);

function LoadProcesses(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
procedure SaveProcesses(filename: String);
function LoadSingleProcess(var F: Textfile):Integer; //returns 0 if succeeded
procedure SaveSingleProcess(var F: Textfile; i: Integer);
procedure SaveSingleProcessRef(var F: TextFile; i: Integer; c: TErgProcess; ref_limitations: Array of TErgLimitation; AlwaysShowDefault: Boolean=true);

function LoadCElements(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
procedure SaveCElements(filename: String);
function LoadSingleCElement(var F: Textfile):Integer; //returns 0 if succeeded
procedure SaveSingleCElement(var F: Textfile; i: Integer);
procedure SaveSingleCElementRef(var F: TextFile; i: Integer; c: TErgCElement; AlwaysShowDefault: Boolean=true);

function LoadModelInfos(filename: String):Integer; //returns 0 if succeeded
procedure SaveModelInfos(filename: String);

procedure WriteConstantsHeader(var F: TextFile);
procedure WriteTracersHeader(var F: TextFile);
procedure WriteAuxiliariesHeader(var F: TextFile);
procedure WriteProcessesHeader(var F: TextFile);
procedure WriteElementsHeader(var F: TextFile);
procedure WriteCElementsHeader(var F: TextFile);
procedure WriteModelinfosHeader(var F: TextFile);

function GenerateIndexes(sort_moving_notmoving: Boolean=false): String;   //fills all the "my..." values in the records
                                    //returns error string or ''
function CheckAuxiliaryOrder: String;
                                    //returns error string or ''
procedure SortAuxiliaries;                                    

procedure FindSourcesSinks(element: Integer; var sources: TStringList; var sinks: TStringList; ConsiderVirtualTracers: Boolean=true);
function ProcessIsConservative(process: Integer; ConsiderVirtualTracers: Boolean=true):Boolean;

function ElementCreatedByProcess_source(element: Integer; process: Integer; ConsiderVirtualTracers: Boolean=true):Real;
function ElementCreatedByProcess_sink(element: Integer; process: Integer; ConsiderVirtualTracers: Boolean=true):Real;
function ElementCreatedByProcess(element: Integer; process: Integer; ConsiderVirtualTracers: Boolean=true):Real;

procedure AutoLimitProcesses;

procedure ApplyLimitations;

procedure GetProcessTypes(var sl: TStringList);
procedure SetDisabledProcessesInactive(sl: TStringList);

function GetInputOutputEquation(c: tErgProcess):String;
function SetInputOutputEquation(s: String; var p: TErgProcess):Boolean;

function SemiItem(var s: String): String;
function StrToIntVertLoc(s:String):Integer;
function IntToStrVertLoc(i:Integer):String;
procedure splitEquation(s: string; var s1: String; var s2: String);
procedure splitSpace(s: string; var s1: String; var s2: String);

procedure initKeywords;

function FinishLoadingAddOn: String;

function GetProcessLimitation(s: String; var prolim: TErgProcessLimitation):Boolean;

procedure AddExtraTagFromText(s1, s2: String; var ExtraTagsRecord: TExtraTagsRecord);
procedure SaveExtraTagsRecord(ExtraTagRecord: TExtraTagsRecord; var F: Textfile; width: Integer = 12);

implementation

uses sysUtils, parser, math;

function SemiItem(var s: String): String;
begin
  if pos(';',s)>0 then
  begin
    result:=copy(s,1,pos(';',s)-1);
    s:=copy(s,pos(';',s)+1,length(s));
  end
  else
  begin
    result:=s;
    s:='';
  end;
end;

function StrToIntVertLoc(s:String):Integer;
begin
  if lowercase(s)='wat' then result:=0
  else if lowercase(s)='sed' then result:=1
  else if lowercase(s)='sur' then result:=2
  else if lowercase(s)='fis' then result:=3
  else result:=StrToInt(s);
end;

function IntToStrVertLoc(i:Integer):String;
begin
  if i=0 then result:='WAT'
  else if i=1 then result:='SED'
  else if i=2 then result:='SUR'
  else if i=3 then result:='FIS'
  else result:=IntToStr(i);
end;

procedure splitEquation(s: string; var s1: String; var s2: String);
var myS: String;
begin
  if pos('!',s)>0 then
    myS:=copy(s,1,pos('!',s)-1)  //! stands for comment
  else
    myS:=s;
  if pos('=',myS)>0 then
  begin
    s1:=copy(myS,1,pos('=',myS)-1);
    s2:=copy(myS,pos('=',myS)+1,length(myS));
  end
  else
  begin
    s1:=myS;
    s2:='';
  end;
  s1:=trim(s1); s2:=trim(s2);
end;

procedure splitSpace(s: string; var s1: String; var s2: String);
var myS: String;
begin
  myS:=s;
  if pos(' ',myS)>0 then
  begin
    s1:=copy(myS,1,pos(' ',myS)-1);
    s2:=copy(myS,pos(' ',myS)+1,length(myS));
  end
  else
  begin
    s1:=myS;
    s2:='';
  end;
  s1:=trim(s1); s2:=trim(s2);
end;

procedure initKeywords;
begin
  setLength(keywords,26); setLength(keywordDescriptions,26);
  keywords[ 0]:='WAT';                       keywordDescriptions[ 0]:='a value for property vertLoc, means "everywhere in the water column" (default)';
  keywords[ 1]:='SED';                       keywordDescriptions[ 1]:='a value for property vertLoc, means "in the sediment only"';
  keywords[ 2]:='SUR';                       keywordDescriptions[ 2]:='a value for property vertLoc, means "at the surface only"';
  keywords[ 3]:='FIS';                       keywordDescriptions[ 3]:='a value for property vertLoc, means "fish-type behaviour" (tracers stored like SED-tracers, but interact with the whole water column)';
  keywords[ 4]:='HARD';                      keywordDescriptions[ 4]:='limitation type (hard), adds a factor "*theta(t-t_min)" to the process turnover';
  keywords[ 5]:='MM';                        keywordDescriptions[ 5]:='limitation type (Michaeles-Menten), adds a factor "*t/(t+t_min)" to the process turnover';
  keywords[ 6]:='MMQ';                       keywordDescriptions[ 6]:='limitation type (quadratic Michaeles-Menten), adds a factor "*t*t/(t*t+t_min*t_min)" to the process turnover';
  keywords[ 7]:='IV';                        keywordDescriptions[ 7]:='limitation type (Ivlev), adds a factor "*(1.0-exp(-t/t_min))" to the process turnover';
  keywords[ 8]:='IVQ';                       keywordDescriptions[ 8]:='limitation type (quadratic Ivlev), adds a factor "*(1.0-exp(-t*t/(t_min*t_min)))" to the process turnover';
  keywords[ 9]:='LIN';                       keywordDescriptions[ 9]:='limitation type (linear), adds a factor "*min(1.0,t/t_min)" to the process turnover';
  keywords[10]:='TANH';                      keywordDescriptions[10]:='limitation type (tangens hyperbolicus), adds a factor "*tanh(t/t_min)" to the process turnover';
  keywords[11]:='cgt_bottomdepth';           keywordDescriptions[11]:='external parameter provided by the host model: Distance between sea surface and bottom of the current cell [m]';
  keywords[12]:='cgt_cellheight';            keywordDescriptions[12]:='external parameter provided by the host model: Height of the current cell [m]';
  keywords[13]:='cgt_current_wave_stress';   keywordDescriptions[13]:='external parameter provided by the host model: Combined bottom stress of current and waves [N/m2]';
  keywords[14]:='cgt_dayofyear';             keywordDescriptions[14]:='external parameter provided by the host model: Julian day in the current year (integer) [days]';
  keywords[15]:='cgt_density';               keywordDescriptions[15]:='external parameter provided by the host model: Density [kg/m3]';
  keywords[16]:='cgt_hour';                  keywordDescriptions[16]:='external parameter provided by the host model: Fractional hours since midnight (0..23.99) [hours]';
  keywords[17]:='cgt_iteration';             keywordDescriptions[17]:='external parameter provided by the host model: Number of current iteration in the iterative loop, starting from 1 (integer) [1]';
  keywords[18]:='cgt_latitude';              keywordDescriptions[18]:='external parameter provided by the host model: Latitude [degrees_north]';
  keywords[19]:='cgt_light';                 keywordDescriptions[19]:='external parameter provided by the host model: Downward light flux (photosynthetically active radiation only) [W/m2]';
  keywords[20]:='cgt_longitude';             keywordDescriptions[20]:='external parameter provided by the host model: Longitude [degrees_east]';
  keywords[21]:='cgt_sali';                  keywordDescriptions[21]:='external parameter provided by the host model: Salinity [psu]';
  keywords[22]:='cgt_temp';                  keywordDescriptions[22]:='external parameter provided by the host model: Potential temperature [°C]';
  keywords[23]:='cgt_timestep';              keywordDescriptions[23]:='external parameter provided by the host model: Time step [days]';
  keywords[24]:='cgt_year';                  keywordDescriptions[24]:='external parameter provided by the host model: Year (integer) [years]';
  keywords[25]:='cgt_diffusivity';           keywordDescriptions[25]:='effective (turbulent) diffusivity in the water column [m2/s]';
end;

function IsStandardChar(c: Char):Boolean;
begin
  if (ord(c)>=ord(char('a'))) and (ord(c)<=ord(char('z'))) then
    result:=true
  else if (ord(c)>=ord(char('A'))) and (ord(c)<=ord(char('Z'))) then
    result:=true
  else if (ord(c)>=ord(char('0'))) and (ord(c)<=ord(char('9'))) then
    result:=true
  else if (ord(c)=ord(char('_'))) or (ord(c)=ord(char('$')))  then
    result:=true
  else
    result:=false;
end;

procedure renameInString(oldName,newName: String; var myString: String);
var
  i: Integer;
  standsAlone: Boolean;
  s: String;
begin
  i:=1;
  while i <= length(myString)-length(oldName)+1 do
  begin
    if lowercase(copy(myString,i,length(oldName))) = lowercase(oldname) then
    begin
      standsAlone:=true;
      if i>1 then
      begin
        s:=copy(myString,i-1,1);
        if IsStandardChar(s[1]) then standsAlone:=false;
      end;
      if i < length(myString)-length(oldName)+1 then
      begin
        s:=copy(myString,i+length(oldName),1);
        if IsStandardChar(s[1]) then standsAlone:=false;
      end;
      if standsAlone then
      begin
        myString:=copy(myString,1,i-1)+newName+copy(myString,i+length(oldName),length(myString));
        i:=i+length(newName)-1;
      end;
    end;
    i:=i+1;
  end;
end;

function containedInString(substr, str: String):Boolean;
var s: String;
begin
  s:=str;
  RenameInString(substr,'',s);
  if s=str then result:=false else result:=true;
end;

procedure AddExtraTagFromText(s1, s2: String; var ExtraTagsRecord: TExtraTagsRecord);
var
  i_tag: Integer;
begin
  if validTagName(s1) then
  begin
    i_tag := ExtraTagsRecord.TagNames.IndexOf(s1) ;
    if i_tag < 0 then //record does not contain tag yet
    begin
      ExtraTagsRecord.TagNames.Add(s1);
      ExtraTagsRecord.TagValues.Add(s2);
    end
    else
      ExtraTagsRecord.TagValues[i_tag] := s2;
  end;
end;


procedure SaveExtraTagsRecord(ExtraTagRecord: TExtraTagsRecord; var F: Textfile; width: Integer = 12);
var
  i_tag, NumberOfSpaces: Integer;
  TextLine: String;
begin
  for i_tag := 0 to ExtraTagRecord.tagNames.Count -1 do
  begin
    if (ExtraTagRecord.tagValues[i_tag] <> '') and
       validTagName(ExtraTagRecord.tagNames[i_tag]) then
    begin
      NumberOfSpaces := Max(0, width - Length(ExtraTagRecord.tagNames[i_tag]));
      TextLine := ExtraTagRecord.tagNames[i_tag] + Space(NumberOfSpaces);
      TextLine := TextLine + '= ' + ExtraTagRecord.tagValues[i_tag];

      WriteLn(F, TextLine);
    end;
  end;
end;


//************************** constants **********************************//

function LoadSingleConstant(var F: TextFile): Integer; //returns 0 if succeeded
var
  s,s1,s2: String;
  i, p: Integer;
begin
  result:=0;
  setLength(Constants,length(Constants)+1);
  i:=length(constants)-1;
  initErgConstant(Constants[length(Constants)-1]);
  s1:=''; s2:='';
  while not EOF(F) do
  begin
      readln(F,s);
      splitEquation(s,s1,s2);
      if s1 <> '' then
      begin
        if (copy(s1,1,1)='*') then            //constant ended
          break
        else if lowercase(s1)='name' then
        begin
          constants[i].name:=s2;
          for p:=0 to i-1 do                  //if constant already exists, update it
            if lowercase(constants[i].name)=lowercase(constants[p].name) then
            begin
              i:=p;
              setLength(constants, length(constants)-1);
            end;
        end
        else if lowercase(s1)='value' then
        begin
          constants[i].valueString:=s2;
          //if pos(';',s2)<=0 then constants[length(constants)-1].value:=StrToFloat(s2);
        end
        else if lowercase(s1)='dependson' then
          constants[i].dependsOn:=lowercase(s2)
        else if lowercase(s1)='minval' then
          constants[i].minval:=StrToFloat(s2)
        else if lowercase(s1)='maxval' then
          constants[i].maxval:=StrToFloat(s2)
        else if lowercase(s1)='varphase' then
          constants[i].varphase:=StrToInt(s2)
        else if lowercase(s1)='description' then
          constants[i].description:=s2
        else if lowercase(s1)='comment' then
          constants[i].comment:=s2
        else if lowercase(s1)='unit' then
          constants[i].valUnit:=s2
        else //undefined variables were set, add as extra tag
          AddExtraTagFromText(s1, s2, constants[i].ExtraTagsRecord);
      end;
  end;    //while ended, now remove possibly empty last entry
  if constants[length(constants)-1].name = '' then
    setLength(constants,length(constants)-1);
end;

function LoadConstants(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  if not loadAddOn then setLength(Constants,0);
  AssignFile(F,filename);
  reset(F);
  result:=0;
    while not EOF(F) do
      result:=max(LoadSingleConstant(F),result);
  closefile(F);
end;

procedure WriteConstantsHeader(var F: TextFile);
begin
  writeln(F,'! BioCGT constants file');
  writeln(F,'! *********************');
  writeln(F,'! properties of constants:');
  writeln(F,'!   name=           code name of constant');
  writeln(F,'!   description=    description including unit, default=""');
  writeln(F,'!   value=          value(s) separated by ";"');
  writeln(F,'!   dependsOn=      whether this constant varies in time and space, default="none", other possible values: "xyz", "xyzt"');
  writeln(F,'!   minval=         minimum value for variation (for inverse modelling)');
  writeln(F,'!   maxval=         maximum value for variation (for inverse modelling)');
  writeln(F,'!   varphase=       phase in which a parameter is allowed to vary (for inverse modelling) (default=1)');
  writeln(F,'!   comment=        comment, e.g. how certain this value is, literature,..., default=""');
  writeln(F,'!   unit=           unit of the constant');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleConstantRef(var F: TextFile; i: Integer; c:TErgConstant; alwaysShowDefault: Boolean=true);
begin
  writeln(F,'name        = '+constants[i].name);
  if (constants[i].value <> c.value) or alwaysShowDefault then
    writeln(F,'value       = '+constants[i].valueString);
  if (constants[i].dependsOn <> c.dependsOn) then
    writeln(F,'dependsOn   = '+constants[i].dependsOn);
  if (constants[i].minval <> c.minval) then
    writeln(F,'minval      = '+FloatToStr(constants[i].minval));
  if (constants[i].maxval <> c.maxval) then
    writeln(F,'maxval      = '+FloatToStr(constants[i].maxval));
  if (constants[i].varphase <> c.varphase) then
    writeln(F,'varphase    = '+IntToStr(constants[i].varphase));
  if (constants[i].description <> c.description) or alwaysShowDefault then
    writeln(F,'description = '+constants[i].description);
  if constants[i].comment <> c.comment then
    writeln(F,'comment     = '+constants[i].comment);
  if constants[i].valUnit <> c.valUnit then
    writeln(F,'unit        = '+constants[i].valUnit);

  //Save extra tags
  SaveExtraTagsRecord(constants[i].ExtraTagsRecord, F);
end;

procedure SaveSingleConstant(var F: TextFile; i: Integer);
var c: TErgConstant;
begin
  InitErgConstant(c);
  SaveSingleConstantRef(F,i,c);
end;


procedure SaveConstants(filename: String);
var F: TextFile;
    i: Integer;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  AssignFile(F,FileName);
  rewrite(F);
    WriteConstantsHeader(F);
    for i:=0 to length(constants)-1 do
    begin
      SaveSingleConstant(F,i);
      writeln(F,'***********************');
    end;
  closefile(F);
end;

//************************** elements **********************************//

function LoadSingleElement(var F: TextFile): Integer; //returns 0 if succeeded
var
  s,s1,s2: String;
  i, p: Integer;
begin
  result:=0;
  setLength(Elements,length(Elements)+1);
  i:=length(elements)-1;
  initErgElement(Elements[length(Elements)-1]);
  s1:=''; s2:='';
  while not EOF(F) do
  begin
      readln(F,s);
      splitEquation(s,s1,s2);
      if s1 <> '' then
      begin
        if (copy(s1,1,1)='*') then            //Element ended
          break
        else if lowercase(s1)='name' then
        begin
          Elements[i].name:=s2;
          for p:=0 to i-1 do                  //if element already exists, update it
            if lowercase(elements[i].name)=lowercase(elements[p].name) then
            begin
              i:=p;
              setLength(elements, length(elements)-1);
              break;
            end;
        end
        else if lowercase(s1)='description' then
          Elements[i].description:=s2
        else if lowercase(s1)='comment' then
          Elements[i].comment:=s2
        else //undefined variables were set, add as extra tag
          AddExtraTagFromText(s1, s2, Elements[i].ExtraTagsRecord);
        Elements[i].isColored:=0;
        Elements[i].hasColoredCopies:=0;
        Elements[i].parentElement:=i;
      end;
  end;    //while ended, now check if a Element with this name exists
  if Elements[length(Elements)-1].name = '' then
    setLength(Elements,length(Elements)-1);
end;

function LoadElements(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  if not loadAddOn then setLength(Elements,0);
  AssignFile(F,filename);
  reset(F);
  result:=0;
    while not EOF(F) do
      result:=max(LoadSingleElement(F),result);
  closefile(F);
end;

procedure WriteElementsHeader(var F: TextFile);
begin
  writeln(F,'! BioCGT elements file');
  writeln(F,'! ********************');
  writeln(F,'! properties of elements:');
  writeln(F,'!   name=           internal name used, e.g. "N"');
  writeln(F,'!   description=    e.g. "nitrogen"');
  writeln(F,'!   comment=        any comments');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleElementRef(var F: TextFile; i: Integer; c:TErgElement; AlwaysShowDefault: Boolean=true);
begin
  writeln(F,'name        = '+elements[i].name);
  if (elements[i].description<>c.description) or AlwaysShowDefault then
    writeln(F,'description = '+elements[i].description);
  if elements[i].comment<>c.comment then
    writeln(F,'comment = '+elements[i].comment);

  //Save extra tags
  SaveExtraTagsRecord(elements[i].ExtraTagsRecord, F);
end;

procedure SaveSingleElement(var F: TextFile; i: Integer);
var c: TErgElement;
begin
  InitErgElement(c);
  SaveSingleElementRef(F,i,c);
end;

procedure SaveElements(filename: String);
var F: TextFile;
    i: Integer;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  AssignFile(F,FileName);
  rewrite(F);
    WriteElementsHeader(F);
    for i:=0 to length(Elements)-1 do
    begin
      SaveSingleElement(F,i);
      writeln(F,'***********************');
    end;
  closefile(F);
end;

//************************** tracers **********************************//

function LoadSingleTracer(var F: TextFile): Integer; //returns 0 if succeeded
var
  s,s1,s2: String;
  i, ii, p: Integer;
begin
  result:=0;
  setLength(Tracers,length(Tracers)+1);
  i:=length(tracers)-1;
  initErgTracer(Tracers[length(Tracers)-1]);
  s1:=''; s2:='';
  while not EOF(F) do
  begin
      readln(F,s);
      splitEquation(s,s1,s2);
      if s1 <> '' then
      begin
        if (copy(s1,1,1)='*') then            //Tracer ended
          break
        else if lowercase(s1)='name' then
        begin
          Tracers[i].name:=s2;
          for p:=0 to i-1 do                  //if tracer already exists, update it
            if lowercase(tracers[i].name)=lowercase(tracers[p].name) then
            begin
              i:=p;
              setLength(tracers, length(tracers)-1);
              break;
            end;
        end
        else if lowercase(s1)='description' then
          Tracers[i].description:=s2
        else if lowercase(s1)='comment' then
          Tracers[i].comment:=s2
        else if lowercase(s1)='solubility' then
          Tracers[i].solubility:=s2
        else if lowercase(s1)='schmidtnumber' then
          Tracers[i].schmidtNumber:=s2
        else if lowercase(s1)='gasname' then
          Tracers[i].gasName:=s2
        else if lowercase(s1)='molarmass' then
          Tracers[i].molarMass:=s2
        else if lowercase(s1)='contents' then
        begin
          SetLength(Tracers[i].contents,StrToInt(s2));
          for ii:=0 to length(Tracers[i].contents)-1 do
          begin
            s1:='';
            while s1='' do
            begin
              readln(F,s);
              splitEquation(s,s1,s2);
            end;
            Tracers[i].contents[ii].element:=s1;
            Tracers[i].contents[ii].amount:=s2;
          end;
        end
        else if lowercase(s1)='verticaldistribution' then
          Tracers[i].verticalDistribution:=s2
        else if lowercase(s1)='vertspeed' then
          Tracers[i].vertSpeed:=s2
        else if lowercase(s1)='vertdiff' then
          Tracers[i].vertDiff:=s2
        else if lowercase(s1)='opacity' then
          Tracers[i].opacity:=s2
        else if (lowercase(s1)='isflat') or (lowercase(s1)='vertloc') then
          Tracers[i].vertLoc:=StrToIntVertLoc(s2)
        else if lowercase(s1)='isinporewater' then
          Tracers[i].isInPorewater:=StrToInt(s2)
        else if lowercase(s1)='ispositive' then
          Tracers[i].isPositive:=StrToInt(s2)
        else if lowercase(s1)='ismixed' then
          Tracers[i].isMixed:=StrToInt(s2)
        else if lowercase(s1)='iscombined' then
          Tracers[i].isCombined:=StrToInt(s2)
        else if lowercase(s1)='isstiff' then
          Tracers[i].isStiff:=StrToInt(s2)
        else if lowercase(s1)='initvalue' then
          Tracers[i].initValue:=s2
        else if lowercase(s1)='useinitvalue' then
          Tracers[i].useInitValue:=StrToInt(s2)
        else if lowercase(s1)='childof' then
          Tracers[i].childOf:=s2
        else if lowercase(s1)='dimension' then
          Tracers[i].dimension:=StrToInt(s2)
        else if lowercase(s1)='masslimits' then
          Tracers[i].massLimits:=s2
        else if lowercase(s1)='atmosdep' then
          Tracers[i].atmosDep:=StrToInt(s2)
        else if lowercase(s1)='riverdep' then
          Tracers[i].riverDep:=StrToInt(s2)
        else if lowercase(s1)='isoutput' then
          Tracers[i].isOutput:=StrToInt(s2)
        else if lowercase(s1)='isactive' then
          Tracers[i].isActive:=StrToInt(s2)
        else if lowercase(s1)='moldiff' then
          Tracers[i].molDiff:=s2
        else if lowercase(s1)='tracerabove' then
          Tracers[i].tracerAbove:=s2
        else if lowercase(s1)='tracerbelow' then
          Tracers[i].tracerBelow:=s2
        else if lowercase(s1)='longname' then
          Tracers[i].longname:=s2
        else if lowercase(s1)='outputunit' then
          Tracers[i].outputunit:=s2
        else if lowercase(s1)='stdname_prefix' then
          Tracers[i].stdname_prefix:=s2
        else if lowercase(s1)='stdname_suffix' then
          Tracers[i].stdname_suffix:=s2
        else //undefined variables were set, add as extra tag
          AddExtraTagFromText(s1, s2, Tracers[i].ExtraTagsRecord);
      end;
  end;  // while ended: check whether longname=''; then set it to name
  if Tracers[i].longname='' then
    Tracers[i].longname:=Tracers[i].name;
  //while ended: now check if a Tracer with this name exists
  if Tracers[length(Tracers)-1].name = '' then
    setLength(Tracers,length(Tracers)-1);
end;

function LoadTracers(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
var F: TextFile; i: integer;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  if loadAddOn=false then
    setLength(Tracers,0);

  AssignFile(F,filename);
  reset(F);
  result:=0;
  while not EOF(F) do
    result := max(LoadSingleTracer(F), result);
  closefile(F);
end;

procedure WriteTracersHeader(var F: TextFile);
begin
  writeln(F,'! BioCGT Tracers file');
  writeln(F,'! *******************');
  writeln(F,'! properties of Tracers:');
  writeln(F,'! name=           variable name in the code');
  writeln(F,'! description=    e.g. "flaggelates"');
  writeln(F,'! atmosDep=       0=no atmospheric depositon of this tracer (default), 1=atmospheric deposition of this tracer');
  writeln(F,'! childOf=        e.g. "flaggelates" for "red_N_in_flaggelates", default="none"');
  writeln(F,'! contents=       number n of elements contained in this tracer, default=0');
  writeln(F,'!                 This line is followed by n lines of this kind:');
  writeln(F,'!   <element> = <quantity>   where <element> is the element name (e.g., "N")');
  writeln(F,'!                            and <quantity> is the quantity of this element in one tracer unit.');
  writeln(F,'!                            Valid for <quantity> are real numbers or names of constants or code expressions, possibly containing auxiliary variables, that can be evaluated without knowing external parameters and tracer concentrations.');
  writeln(F,'! dimension=      how many instances of this tracer exist, e.g. a tracer named "cod" with dimension=2 exists as "cod_$cod" which is "cod_1" and "cod_2", default=0');
  writeln(F,'! gasName=        name of an auxiliary variable containing the concentration [mol/kg] of the dissolved gas in the surface cell, e.g. "co2" for tracer "dic". Default="" meaning it is the tracer concentration itself');
  writeln(F,'! initValue=      initial value, default="0.0", set "useInitValue" to 1 to use it');
  writeln(F,'! isActive=       1=active (default); 0=virtual tracer to check element conservation');
  writeln(F,'! isCombined=     1=combined tracer that accumulates several virtual tracers (isActive=false) in one variable, its contents are tracers rather than elements; 0=not (default)');
  writeln(F,'! isMixed=        1=mix with neighbour cells if negative, 0=do not, default=0, only applicable to tracers with vertLoc=WAT');
  writeln(F,'! isOutput=       1=occurs as output in model results (default); 0=only internal use');
  writeln(F,'! isInPorewater=  1=dissolved species which is also present in porewater, 0=not (particulate species) (default)');  
  writeln(F,'! isPositive=     0=may be negative, 1=always positive (default)');
  writeln(F,'! isStiff=        0=not stiff (default); 1=stiff, use Patankar method if concentration declines; 2=stiff, always use modified Patankar method');
  writeln(F,'! massLimits=     semicolon-seperated string of (dimension-1) FORTRAN expressions for mass limits for stage-resolving models [mmol], default=""');
  writeln(F,'! molDiff=        molecular diffusivity in pore water [m2/s], use the name of a vertLoc=SED auxiliary variable, default="0.0"');
  writeln(F,'! opacity=        fortran formula for opacity [m2/mol] (for light limitation), default="0"');
  writeln(F,'! riverDep=       0=no river depositon of this tracer, 1=river deposition of this tracer (default)');
  writeln(F,'! schmidtNumber=  name of an auxiliary variable describing the Schmidt number [1] for gasses which flow through the surface, default="0"');
  writeln(F,'! solubility=     name of an auxiliary variable describing the solubility [mol/kg/Pa] for gasses which flow through the surface, default="0"');
  writeln(F,'! useInitValue=   1=use initValue as initial concentration, 0=do not (load initial concentration from file) (default)');
  writeln(F,'! verticalDistribution= Name of an auxiliary variable proportional to which the vertical distribution of the tracer is assumed. Relevant for vertLoc=FIS only. Default="1.0"');
  writeln(F,'! vertDiff=       fortran formula for vertical diffusivity [m2/s], default="0"');
  writeln(F,'! vertLoc=        WAT=everywhere in the water column (default), SED=in sediment only, SUR=at surface only, FIS=fish-type behaviour');
  writeln(F,'! vertSpeed=      fortran formula for vertical speed [m/day], default="0"');
  writeln(F,'! tracerAbove=    name of corresponding tracer in the compartment above, default="none" (e.g. "det" (in water column) if we are "benth_det")');
  writeln(F,'! tracerBelow=    name of corresponding tracer in the compartment below, default="none" (e.g. "nit" (in water column) if we are "surf_nit")');
  writeln(F,'! longname=       longname of the tracer (i.e. nitrate for a tracer named nit); used e.g. for the longname attributes in netCDF files');
  writeln(F,'! outputUnit=     unit of the tracer, which should be printed by the <unit> tag; e.g. for the unit attributes in netCDF files');
  writeln(F,'! stdname_prefix= prefix for a CF-convention standardname (inserted by the <standardname> tag), which is build by stdname_prefix+longname+stdname_suffix');
  writeln(F,'! stdname_suffix= suffix for a CF-convention standardname (inserted by the <standardname> tag), which is build by stdname_prefix+longname+stdname_suffix');
  writeln(F,'! comment=        e.g. "represents a certain kind of phytoplankton", default=""');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleTracerRef(var F: TextFile; i: Integer; c: TErgTracer; AlwaysShowDefault: Boolean=true);
var
  j: Integer;
  writeContents: Boolean;
//  writeExtraTags: Boolean;
begin
  writeln(F,'name              = '+Tracers[i].name);
  if (Tracers[i].description <> c.description) or AlwaysShowDefault then
    writeln(F,'description       = '+Tracers[i].description);
  if Tracers[i].comment <> c.comment then
    writeln(F,'comment           = '+Tracers[i].comment);

  writeContents := false;
  if length(tracers[i].contents) <> length(c.contents) then
    writeContents := true
  else
  begin
    for j:=0 to length(tracers[i].contents)-1 do
    begin
      if tracers[i].contents[j].element <> c.contents[j].element then
        writeContents := true;
      if tracers[i].contents[j].amount <> c.contents[j].amount then
        writeContents := true;
    end;
  end;
  if WriteContents then
  begin
    writeln(F,'contents          = '+IntToStr(length(Tracers[i].contents)));
    for j:=0 to length(Tracers[i].contents)-1 do
    begin
      writeln(F,'  '+Tracers[i].contents[j].element+' = '+ Tracers[i].contents[j].amount);
    end;
  end;

  if Tracers[i].isOutput <> c.isOutput then
    writeln(F,'isOutput          = '+IntToStr(Tracers[i].isOutput));
  if Tracers[i].verticalDistribution <> c.verticalDistribution then
    writeln(F,'verticalDistribution='+Tracers[i].verticalDistribution);
  if Tracers[i].vertSpeed <> c.vertSpeed then
    writeln(F,'vertSpeed         = '+Tracers[i].vertSpeed);
  if Tracers[i].vertDiff <> c.vertDiff then
    writeln(F,'vertDiff          = '+Tracers[i].vertDiff);
  if Tracers[i].opacity <> c.opacity then
    writeln(F,'opacity           = '+Tracers[i].opacity);
  if Tracers[i].vertLoc <> c.vertLoc then
    writeln(F,'vertLoc           = '+IntToStrVertLoc(Tracers[i].vertLoc));
  if Tracers[i].isInPorewater <> c.isInPorewater then
    writeln(F,'isInPorewater     = '+IntToStr(Tracers[i].isInPorewater));
  if Tracers[i].isPositive <> c.isPositive then
    writeln(F,'isPositive        = '+IntToStr(Tracers[i].isPositive));
  if Tracers[i].isMixed <> c.isMixed then
    writeln(F,'isMixed           = '+IntToStr(Tracers[i].isMixed));
  if Tracers[i].isCombined <> c.isCombined then
    writeln(F,'isCombined        = '+IntToStr(Tracers[i].isCombined));
  if Tracers[i].isStiff <> c.isStiff then
    writeln(F,'isStiff           = '+IntToStr(Tracers[i].isStiff));
  if Tracers[i].isActive <> c.isActive then
    writeln(F,'isActive          = '+IntToStr(Tracers[i].isActive));
  if Tracers[i].childOf <> c.childOf then
    writeln(F,'childOf           = '+Tracers[i].childOf);
  if Tracers[i].initValue <> c.initValue then
    writeln(F,'initValue         = '+Tracers[i].initValue);
  if Tracers[i].useInitValue <> c.useInitValue then
    writeln(F,'useInitValue      = '+IntToStr(Tracers[i].useInitValue));
  if Tracers[i].isOutput <> c.isOutput then
    writeln(F,'isOutput          = '+IntToStr(Tracers[i].isOutput));
  if Tracers[i].solubility <> c.solubility then
    writeln(F,'solubility        = '+Tracers[i].solubility);
  if Tracers[i].schmidtNumber <> c.schmidtNumber then
    writeln(F,'schmidtNumber     = '+Tracers[i].schmidtNumber);
  if Tracers[i].gasName <> c.gasName then
    writeln(F,'gasName           = '+Tracers[i].gasName);
  if Tracers[i].molarMass <> c.molarMass then
    writeln(F,'molarMass         = '+Tracers[i].molarMass);
  if Tracers[i].dimension <> c.dimension then
    writeln(F,'dimension         = '+IntToStr(Tracers[i].dimension));
  if Tracers[i].massLimits <> c.massLimits then
    writeln(F,'massLimits        = '+Tracers[i].massLimits);
  if Tracers[i].atmosDep <> c.atmosDep then
    writeln(F,'atmosDep          = '+IntToStr(Tracers[i].atmosDep));
  if Tracers[i].riverDep <> c.riverDep then
    writeln(F,'riverDep          = '+IntToStr(Tracers[i].riverDep));
  if Tracers[i].molDiff <> c.molDiff then
    writeln(F,'molDiff           = '+Tracers[i].molDiff);
  if Tracers[i].tracerAbove <> c.tracerAbove then
    writeln(F,'tracerAbove       = '+Tracers[i].tracerAbove);
  if Tracers[i].tracerBelow <> c.tracerBelow then
    writeln(F,'tracerBelow       = '+Tracers[i].tracerBelow);
  if Tracers[i].longname <> Tracers[i].name then  // note: different comparison!!!
    writeln(F,'longname          = '+Tracers[i].longname);
  if Tracers[i].outputUnit <> c.outputUnit then
    writeln(F,'outputUnit        = '+Tracers[i].outputUnit);
  if Tracers[i].stdname_prefix <> c.stdname_prefix then
    writeln(F,'stdname_prefix    = '+Tracers[i].stdname_prefix);
  if Tracers[i].stdname_suffix <> c.stdname_suffix then
    writeln(F,'stdname_suffix    = '+Tracers[i].stdname_suffix);

  //Save extra tags
  SaveExtraTagsRecord(Tracers[i].ExtraTagsRecord, F, 18);
end;

procedure SaveSingleTracer(var F: TextFile; i: Integer);
var
  c: TErgTracer;
begin
  InitErgTracer(c);
  //InitErgExtraTags(c.extraTags); //sets number of tags to 0
  SaveSingleTracerRef(F,i,c);
end;

procedure SaveTracers(filename: String);
var F: TextFile;
    i: Integer;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  AssignFile(F,FileName);
  rewrite(F);
    WriteTracersHeader(F);
    for i:=0 to length(Tracers)-1 do
    begin
      SaveSingleTracer(F,i);
      writeln(F,'***********************');
    end;
  closefile(F);
end;

//************************** auxiliaries **********************************//

function LoadSingleAuxiliary(var F: TextFile): Integer; //returns 0 if succeeded
var
  s,s1,s2: String;
  i, p: Integer;
begin
  result:=0;
  setLength(Auxiliaries,length(Auxiliaries)+1);
  i:=length(auxiliaries)-1;
  initErgAuxiliary(Auxiliaries[length(Auxiliaries)-1]);
  s1:=''; s2:='';
  while not EOF(F) do
  begin
      readln(F,s);
      splitEquation(s,s1,s2);
      if s1 <> '' then
      begin
        if (copy(s1,1,1)='*') then            //Auxiliary ended
          break
        else if lowercase(s1)='name' then
        begin
          auxiliaries[i].name:=s2;
          for p:=0 to i-1 do                  //if auxiliary already exists, update it
            if lowercase(auxiliaries[i].name)=lowercase(auxiliaries[p].name) then
            begin
              i:=p;
              setLength(auxiliaries, length(auxiliaries)-1);
            end;
        end
        else if lowercase(s1)='temp1' then
          auxiliaries[i].temp[1]:=s2
        else if lowercase(s1)='temp2' then
          auxiliaries[i].temp[2]:=s2
        else if lowercase(s1)='temp3' then
          auxiliaries[i].temp[3]:=s2
        else if lowercase(s1)='temp4' then
          auxiliaries[i].temp[4]:=s2
        else if lowercase(s1)='temp5' then
          auxiliaries[i].temp[5]:=s2
        else if lowercase(s1)='temp6' then
          auxiliaries[i].temp[6]:=s2
        else if lowercase(s1)='temp7' then
          auxiliaries[i].temp[7]:=s2
        else if lowercase(s1)='temp8' then
          auxiliaries[i].temp[8]:=s2
        else if lowercase(s1)='temp9' then
          auxiliaries[i].temp[9]:=s2
        else if lowercase(s1)='formula' then
          auxiliaries[i].formula:=s2
        else if lowercase(s1)='calcafterprocesses' then
          auxiliaries[i].calcAfterProcesses:=strToInt(s2)
        else if lowercase(s1)='iterations' then
          auxiliaries[i].iterations:=StrToInt(s2)
        else if lowercase(s1)='iterinit' then
          auxiliaries[i].iterInit:=s2
        else if lowercase(s1)='isusedelsewhere' then
          auxiliaries[i].isUsedElsewhere:=strToInt(s2)
        else if lowercase(s1)='iszgradient' then
          auxiliaries[i].isZGradient:=StrToInt(s2)
        else if lowercase(s1)='iszintegral' then
          auxiliaries[i].isZIntegral:=StrToInt(s2)
        else if lowercase(s1)='description' then
          auxiliaries[i].description:=s2
        else if lowercase(s1)='comment' then
          auxiliaries[i].comment:=s2
        else if lowercase(s1)='unit' then
          auxiliaries[i].valUnit:=s2
        else if (lowercase(s1)='isflat') or (lowercase(s1)='vertloc') then
          auxiliaries[i].vertLoc:=strToIntVertLoc(s2)
        else if lowercase(s1)='isoutput' then
          auxiliaries[i].isOutput:=strToInt(s2)
        else //undefined variables were set, add as extra tag
          AddExtraTagFromText(s1, s2, auxiliaries[i].ExtraTagsRecord);
      end;
  end;    //while ended, now check if a Auxiliary with this name exists
  if Auxiliaries[length(Auxiliaries)-1].name = '' then
    setLength(Auxiliaries,length(Auxiliaries)-1);
end;

function LoadAuxiliaries(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  if not loadAddOn then setLength(Auxiliaries,0);
  AssignFile(F,filename);
  reset(F);
  result:=0;
    while not EOF(F) do
      result:=max(LoadSingleAuxiliary(F),result);
  closefile(F);
end;

procedure WriteauxiliariesHeader(var F: TextFile);
begin
  writeln(F,'! BioCGT auxiliaries file');
  writeln(F,'! ***********************');
  writeln(F,'! properties of auxiliaries:  (auxiliary values that are calculated)');
  writeln(F,'! name=               variable name in the code');
  writeln(F,'! description=        e.g. "absolute temperature [K]" default=""');
  writeln(F,'! temp1= ... temp9=   for calculating a temporary value which appears in the formula. e.g. temp1=no3*no3 temp2=no3limit*no3limit formula=temp1/(temp1+temp2), default=""');
  writeln(F,'! formula=            formula as it appears in the code');
  writeln(F,'! calcAfterProcesses= 1=calculate this auxiliary after all process rates are known, default=0');
  writeln(F,'! iterations=         how often this auxiliary variable is calculated in an iterative loop, default=0');
  writeln(F,'! iterInit=           the initial value in the iterative loop, default="0.0"');
  writeln(F,'! isOutput=           1=occurs as output in model results; 0=internal use only (default)');
  writeln(F,'! isUsedElsewhere=    1=make the value of this auxiliary accessible from outside the biological model (e.g. use a "diagnostic tracer" in MOM4); 0=internal use only (default)');
  writeln(F,'! isZGradient=        1=is a vertical gradient of a tracer, 0=is not (default). If 1, "formula" must be the name of the tracer, which must have vertLoc=WAT. isZGradient=1 requires vertLoc=WAT.');
  writeln(F,'! isZIntegral=        1=is a vertical integral (of value times density) of a tracer or an auxiliary variable, 0=is not (default). If 1 "formula" must be the name of the tracer, which must have vertLoc=WAT. isZIntegral=1 requires vertLoc=SED.');
  writeln(F,'! vertLoc=            WAT=z-dependent (default), SED=in the bottom cell only, SUR=in the surface cell only');
  writeln(F,'! comment=            e.g. a literature reference, default=""');
  writeln(F,'! unit=               unit of the auxiliary variable');
  writeln(F,'!');
  writeln(F,'! All entries with the same value of calcAfterProcesses are calculated in given order.');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleAuxiliaryRef(var F: TextFile; i: Integer; c: TErgAuxiliary; AlwaysShowDefault: Boolean=true);
var
    j: Integer;
begin
  writeln(F,  'name               = '+auxiliaries[i].name);
  for j:=1 to 9 do
    if auxiliaries[i].temp[j] <> c.temp[j] then
      writeln(F,'temp'+IntToStr(j)+'              = '+auxiliaries[i].temp[j]);
  if (auxiliaries[i].formula <> c.formula) or AlwaysShowDefault then
    writeln(F,  'formula            = '+auxiliaries[i].formula);
  if auxiliaries[i].calcAfterProcesses <> c.calcAfterProcesses then
    writeln(F,'calcAfterProcesses = '+IntToStr(auxiliaries[i].calcAfterProcesses));
  if auxiliaries[i].iterations <> c.iterations then
    writeln(F,'iterations         = '+IntToStr(auxiliaries[i].iterations));
  if auxiliaries[i].iterInit <> c.iterInit then
    writeln(F,'iterInit           = '+auxiliaries[i].iterInit);
  if (auxiliaries[i].description <> c.description) or AlwaysShowDefault then
    writeln(F,'description        = '+auxiliaries[i].description);
  if auxiliaries[i].comment <> c.comment then
    writeln(F,'comment            = '+auxiliaries[i].comment);
  if auxiliaries[i].valUnit <> c.valUnit then
    writeln(F,'unit               = '+auxiliaries[i].valUnit);
  if auxiliaries[i].vertLoc <> c.vertLoc then
    writeln(F,'vertLoc            = '+IntToStrVertLoc(auxiliaries[i].vertLoc));
  if auxiliaries[i].isOutput <> c.isOutput then
    writeln(F,'isOutput           = '+IntToStr(auxiliaries[i].isOutput));
  if auxiliaries[i].isUsedElsewhere <> c.isUsedElsewhere then
    writeln(F,'isUsedElsewhere    = '+IntToStr(auxiliaries[i].isUsedElsewhere));
  if auxiliaries[i].isZGradient <> c.isZGradient then
    writeln(F,'isZGradient        = '+IntToStr(auxiliaries[i].isZGradient));
  if auxiliaries[i].isZIntegral <> c.isZIntegral then
    writeln(F,'isZIntegral        = '+IntToStr(auxiliaries[i].isZIntegral));

  //Save extra tags
  SaveExtraTagsRecord(auxiliaries[i].ExtraTagsRecord, F, 19);
end;

procedure SaveSingleAuxiliary(var F: TextFile; i: Integer);
var
  c: TErgAuxiliary;
begin
  InitErgAuxiliary(c);
  SaveSingleAuxiliaryRef(F,i,c);
end;

procedure SaveAuxiliaries(filename: String);
var F: TextFile;
    i: Integer;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  AssignFile(F,FileName);
  rewrite(F);
    WriteAuxiliariesHeader(F);
    for i:=0 to length(Auxiliaries)-1 do
    begin
      SaveSingleAuxiliary(F,i);
      writeln(F,'***********************');
    end;
  closefile(F);
end;

//************************** Processes **********************************//

function SetInputOutputEquation(s: String; var p: TErgProcess):Boolean;
// Interprets a process equation that is of the following shape:
// value1 * tracer1 + value2 * tracer2 + ... -> value3 * tracer3 + ...
// value1, value2, ... are expressions that the parser can interpret,
//   in the most simple case they are constants from constants.txt or real constants.
// tracer1, tracer2, ... are names of tracers from tracers.txt
// The corresponding arrays p.input and p.output are filled accordingly.
// The function returns true if succeeded and false otherwise.
var
  s1, s2, snew, scurrent, symbol, svalue, stracer: String;
begin
  result:=true;
  if pos('->',s)<=0 then
    result:=false
  else
  begin
    s1:=trim(copy(s,1,pos('->',s)-1));          //left  hand side of the equation = input
    s2:=trim(copy(s,pos('->',s)+2,length(s)));  //right hand side of the equation = output

    //first, inputs
    setLength(p.input,0);
    while FindLastPlusMinus(s1,snew,scurrent,symbol) and (result=true) do
    begin
      //get the last expression valueN * tracerN
      if symbol <> '+' then //only + are allowed
        result:=false
      else
      begin
        if FindLastTimesOver(scurrent,svalue,stracer,symbol) then
        begin
          if symbol <> '*' then
            result:=false;
        end
        else  //no '*' is found, so we assume that "valueN" is 1
        begin
          sValue:='1.0';
          sTracer:=sCurrent;
        end;
        if result=true then //everything is fine until now.
        begin
          setLength(p.input,length(p.input)+1);
          p.input[length(p.input)-1].tracer:=trim(sTracer);
          p.input[length(p.input)-1].amount:=trim(sValue);
        end;
      end;
      s1:=snew; //remove the last entry from the left-hand side of the equation
    end;
    if (s1 <> '') and (result=true) then //get the first entry
    begin
      if FindLastTimesOver(s1,svalue,stracer,symbol) then
      begin
        if symbol <> '*' then
          result:=false;
      end
      else  //no '*' is found, so we assume that "valueN" is 1
      begin
        sValue:='1.0';
        sTracer:=s1;
      end;
      if result=true then //everything is fine until now.
      begin
        setLength(p.input,length(p.input)+1);
        p.input[length(p.input)-1].tracer:=trim(sTracer);
        p.input[length(p.input)-1].amount:=trim(sValue);
      end;
    end;

    //second, outputs, same procedure
    setLength(p.output,0);
    while FindLastPlusMinus(s2,snew,scurrent,symbol) and (result=true) do
    begin
      //get the last expression valueN * tracerN
      if symbol <> '+' then //only + are allowed
        result:=false
      else
      begin
        if FindLastTimesOver(scurrent,svalue,stracer,symbol) then
        begin
          if symbol <> '*' then
            result:=false;
        end
        else  //no '*' is found, so we assume that "valueN" is 1
        begin
          sValue:='1.0';
          sTracer:=sCurrent;
        end;
        if result=true then //everything is fine until now.
        begin
          setLength(p.output,length(p.output)+1);
          p.output[length(p.output)-1].tracer:=trim(sTracer);
          p.output[length(p.output)-1].amount:=trim(sValue);
        end;
      end;
      s2:=snew; //remove the last entry from the left-hand side of the equation
    end;
    if (s2 <> '') and (result=true) then //get the first entry
    begin
      if FindLastTimesOver(s2,svalue,stracer,symbol) then
      begin
        if symbol <> '*' then
          result:=false;
      end
      else  //no '*' is found, so we assume that "valueN" is 1
      begin
        sValue:='1.0';
        sTracer:=s2;
      end;
      if result=true then //everything is fine until now.
      begin
        setLength(p.output,length(p.output)+1);
        p.output[length(p.output)-1].tracer:=trim(sTracer);
        p.output[length(p.output)-1].amount:=trim(sValue);
      end;
    end;
  end;
end;

function GetInputOutputEquation(c: TErgProcess):String;
var
  s: String;
  j: Integer;
begin
      s:='';
      for j:=0 to length(c.input)-1 do
      begin
        if j>0 then s:=s + ' + ';
        if (c.input[j].amount='1') or (c.input[j].amount='1.0') then
          s:=s+c.input[j].tracer
        else
          s:=s+c.input[j].amount+'*'+c.input[j].tracer;
      end;
      s:=s+' -> ';
      for j:=0 to length(c.output)-1 do
      begin
        if j>0 then s:=s + ' + ';
        if (c.output[j].amount='1') or (c.output[j].amount='1.0') then
          s:=s+c.output[j].tracer
        else
          s:=s+c.output[j].amount+'*'+c.output[j].tracer;
      end;
      result:=s;
end;

function GetProcessLimitation(s: String; var prolim: TErgProcessLimitation):Boolean;
var
  limitationType: String;
  tracer: String;
  t, tracerNum: Integer;
  value: String;
  elseProcess: String;
begin
  //first, seek if tracer < value is specified (inverse limitation)
  if pos('<',s)>0 then
  begin
    s:=StringReplace(s,'<','>',[rfReplaceAll,rfIgnoreCase]);
    prolim.tracerIsSmall := 1;
  end
  else
    prolim.tracerIsSmall := 0;
  //find limitation type (HARD, MM, IV, ...)
  s:=trim(s);
  limitationType := uppercase(copy(s,1,pos(' ',s)-1));
  s := trim(copy(s,pos(' ',s)+1,length(s)));
  //find tracer name
  tracerNum:=-1;
  tracer := trim(lowercase(copy(s,1,pos('>',s)-1)));
  s := trim(copy(s,pos('>',s)+1,length(s)));
  for t:=0 to length(tracers)-1 do
    if lowercase(tracers[t].name)=tracer then
      tracerNum:=t;
  if TracerNum = -1 then
    result:=false //tracer not found
  else
  begin
    if pos(' else ',lowercase(s))>0 then
    begin
      value := trim(copy(s,1,pos(' else ',lowercase(s))-1));
      elseProcess := trim(copy(s,pos(' else ',lowercase(s))+6, length(s)));
    end
    else
    begin
      value:=s;
      elseProcess:='';
    end;
    //search if this limitation function already exists
    prolim.limitationNum:=-1;
    for t:=0 to length(limitations)-1 do
      if (limitations[t].limitationType=limitationType)
         and (limitations[t].tracer=tracer)
         and (lowercase(limitations[t].value)=lowercase(value)) then
        prolim.limitationNum:=t;
    if prolim.limitationNum=-1 then //this limitation function did not yet exist, so create it
    begin
      setLength(limitations,length(limitations)+1);
      limitations[length(limitations)-1].limitationType:=limitationType;
      limitations[length(limitations)-1].tracer:=tracer;
      limitations[length(limitations)-1].value:=value;
      prolim.limitationNum:=length(limitations)-1;
    end;
    prolim.elseProcess:=elseProcess;
    result:=true;
  end;
end;

function LoadSingleProcess(var F: TextFile):Integer;
// Loads a single process from the open file F, until a line starts with *
// If its name already exists, it is only modified, else it is appended to the process list.
  //returns 0 if succeeded
  //        2 if undefined variables were set
  //        3 if an error in the process equation occurred
  //        4 if an error occurred in the limitations
var
  s,s1,s2: String;
  i,ii,p: Integer;
  prolim:TErgProcessLimitation;
  IsFirstLimitation: Boolean;
begin
  result:=0;
  setLength(Processes,length(processes)+1);
  i:=length(processes)-1;
  initErgProcess(Processes[length(Processes)-1]);
  s1:=''; s2:='';
  isFirstLimitation:=true;
  while not EOF(F) do
  begin
    readln(F,s);
    splitEquation(s,s1,s2);
    if s1 <> '' then
    begin
      if (copy(s1,1,1)='*') then       //Process ended
        break
      else if lowercase(s1)='name' then
      begin
        Processes[i].name:=s2;
          for p:=0 to i-1 do                  //if process already exists, update it
            if lowercase(processes[i].name)=lowercase(processes[p].name) then
            begin
              i:=p;
              setLength(processes, length(processes)-1);
            end;
        end
      else if lowercase(s1)='description' then
        Processes[i].description:=s2
      else if lowercase(s1)='comment' then
        Processes[i].comment:=s2
      else if (lowercase(s1)='rate') or (lowercase(s1)='turnover') then
        Processes[i].turnover:=s2
      else if lowercase(s1)='comment' then
        Processes[i].comment:=s2
      else if lowercase(s1)='input' then
      begin
        SetLength(Processes[i].input,StrToInt(s2));
        for ii:=0 to length(Processes[i].input)-1 do
        begin
          s1:='';
          while s1='' do
          begin
            readln(F,s);
            splitEquation(s,s1,s2);
          end;
          Processes[i].input[ii].tracer:=s1;
          Processes[i].input[ii].amount:=s2;
        end;
      end
      else if lowercase(s1)='output' then
      begin
        SetLength(Processes[i].output,StrToInt(s2));
        for ii:=0 to length(Processes[i].output)-1 do
        begin
          s1:='';
          while s1='' do
          begin
            readln(F,s);
            splitEquation(s,s1,s2);
          end;
          Processes[i].output[ii].tracer:=s1;
          Processes[i].output[ii].amount:=s2;
        end;
      end
      else if lowercase(s1)='equation' then
      begin
        if SetInputOutputEquation(s2,Processes[i])=false then
          result:=3;
      end
      else if lowercase(s1)='repaint' then
      begin
        SetLength(Processes[i].repaint,StrToInt(s2));
        for ii:=0 to length(Processes[i].repaint)-1 do
        begin
          s1:='';
          while s1='' do
          begin
            readln(F,s);
            splitEquation(s,s1,s2);
          end;
          splitSpace(s1,s,s1);
          Processes[i].repaint[ii].oldColor:=s;
          Processes[i].repaint[ii].element:=s1;
          Processes[i].repaint[ii].newColor:=s2;
        end;
      end
      else if lowercase(s1)='limitation' then
      begin
        if IsFirstLimitation then
        begin
          setLength(Processes[i].limitations,0);
          IsFirstLimitation:=false;
        end;
        if trim(lowercase(copy(s,pos('=',s)+1,length(s))))='none' then //in an add-on, all limitations were removed
        begin
          setLength(Processes[i].limitations,0);
        end
        else if GetProcessLimitation(s2,prolim) then
        begin
          setLength(Processes[i].limitations,length(Processes[i].limitations)+1);
          Processes[i].limitations[length(Processes[i].limitations)-1]:=prolim;
        end
        else
          result:=4;
      end
      else if lowercase(s1)='feedingefficiency' then
        Processes[i].feedingEfficiency:=s2
      else if (lowercase(s1)='isflat') or (lowercase(s1)='vertloc') then
        Processes[i].vertLoc:=StrToIntVertLoc(s2)
      else if lowercase(s1)='processtype' then
        Processes[i].processType:=s2
      else if lowercase(s1)='isoutput' then
        Processes[i].isOutput:=StrToInt(s2)
      else if lowercase(s1)='isactive' then
        Processes[i].isActive:=StrToInt(s2)
      else //undefined variables were set, add as extra tag
      begin
        AddExtraTagFromText(s1, s2, Processes[i].ExtraTagsRecord);
      end;
    end;
  end;  //while ended, now check if a process with this name exists
  if Processes[length(Processes)-1].name = '' then
    setLength(processes,length(processes)-1);
end;

function LoadProcesses(filename: String; loadAddOn: Boolean=false):Integer;
  //returns 0 if succeeded
  //        1 if file could not be loaded
  //        2 if undefined variables were set
  //        3 if an error in the process equation occurred
  //        4 if an error occurred in the limitations
var F: TextFile;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  if not loadAddOn then setLength(limitations,0);
  if not loadAddOn then setLength(Processes,0);
  AssignFile(F,filename);
  reset(F);
  result:=0;
    while not EOF(F) do
      result:=max(LoadSingleProcess(F),result);
  closefile(F);
end;

procedure WriteProcessesHeader(var F: TextFile);
begin
  writeln(F,'! BioCGT processes file');
  writeln(F,'! *********************');
  writeln(F,'! properties of processes:');
  writeln(F,'! name=           fortran variable name for the rate');
  writeln(F,'! description=    e.g. "grazing of zooplankton"');
  writeln(F,'! turnover=       fortran formula for calculating the process turnover [mol/kg or mol/m2]');
  writeln(F,'! equation=       equation which, like a chemical equation, lists reaction agents and products of this process.');
  writeln(F,'!                   example: t_no3 + 1/16*t_po4 -> t_lpp');
  writeln(F,'!                   tracers to the left of the "->" are consumed, tracers to the right of the "->" are produced by this process.');
  writeln(F,'! feedingEfficiency= name of an auxiliary variable (values 0..1) which tells how much of the food in a certain depth is accessible for the predator with vertLoc=FIS. Relevant for vertLoc=FIS only. Default="1.0"');  
  writeln(F,'! isActive=       1=active (default); 0=process is switched off');
  writeln(F,'! isOutput=       1=occurs as output in model results; 0=internal use only (default)');
  writeln(F,'! limitation=     TYPE tracer > value else otherProcess');
  writeln(F,'! limitation=     TYPE tracer < value else otherProcess');
  writeln(F,'!                   TYPE = HARD (theta function), MM (Michaelis-Menten), MMQ (quadratic Michaelis-Menten), IV (Ivlev), IVQ (quadratic Ivlev), LIN (linear), TANH (tangens hyperbolicus)');
  writeln(F,'!                   tracer = name of tracer that needs to be present');
  writeln(F,'!                   value = value that needs to be exceeded, may also be a constant or auxiliary');
  writeln(F,'!                   otherProcess = process that takes place instead if this process gets hampered by the availability of "tracer"');
  writeln(F,'!                 several of these lines may exist, the order of them may be relevant for the rate of "otherProcess".');
  writeln(F,'! processType=    type of process, e.g. "propagation", default="standard"');
  writeln(F,'! repaint=        number n of repainting actions to be done by the process, default=0');
  writeln(F,'!                 This line is followed by n lines of this kind:');
  writeln(F,'!   <oldColor> <element> = <newColor>    e.g.: "all  N   = blue "');
  writeln(F,'!                                              "blue P   = none "');
  writeln(F,'!                                              "red  all = green"');
  writeln(F,'! vertLoc=        WAT=z-dependent (default), SED=in the sediment only, SUR=in the surface only, FIS=fish-type behaviour');
  writeln(F,'! comment=        comment, default=""');
  writeln(F,'!');
  writeln(F,'! Process rates are calculated in the given order.');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleProcessRef(var F: TextFile; i: Integer; c: TErgProcess; ref_limitations: Array of TErgLimitation; AlwaysShowDefault: Boolean=true);
var
  j, k: Integer;
  s, s1: String;
  differs, found: Boolean;
begin
  writeln(F,'name        = '+Processes[i].name);
  if (Processes[i].description <> c.description) or AlwaysShowDefault then
    writeln(F,'description = '+Processes[i].description);
  if (Processes[i].turnover <> c.turnover) or AlwaysShowDefault then
    writeln(F,'turnover    = '+Processes[i].turnover);
  if Processes[i].comment <> c.comment then
    writeln(F,'comment     = '+Processes[i].comment);

  //now, create equation string
  s:=GetInputOutputEquation(Processes[i]);
  //check if input and output are the same, but they do not need to be in the same order
  if (length(processes[i].input)<>length(c.input)) or (length(processes[i].output)<>length(c.output)) then
    differs:=true
  else
  begin
    differs:=false;
    for j:=0 to length(processes[i].input)-1 do
    begin
      found:=false;
      for k:=0 to length(processes[i].input)-1 do
        if (lowercase(processes[i].input[j].tracer)=lowercase(c.input[k].tracer)) and
           (lowercase(processes[i].input[j].amount)=lowercase(c.input[k].amount)) then
          found:=true;
      if not found then differs:=true;
    end;
    if not differs then
      for j:=0 to length(processes[i].output)-1 do
      begin
        found:=false;
        for k:=0 to length(processes[i].output)-1 do
          if (lowercase(processes[i].output[j].tracer)=lowercase(c.output[k].tracer)) and
             (lowercase(processes[i].output[j].amount)=lowercase(c.output[k].amount)) then
            found:=true;
        if not found then differs:=true;
      end;
  end;
  if differs or AlwaysShowDefault then
    writeln(F,'equation    = '+s);
  //check if repaint block needs to be written
  if length(Processes[i].repaint) <> length(c.repaint) then differs:=true
  else
  begin
    differs:=false;
    for j:=0 to length(Processes[i].repaint)-1 do
    begin
      s:='  '+Processes[i].repaint[j].oldColor+' '+Processes[i].repaint[j].element+' = '+ Processes[i].repaint[j].newColor;
      s1:='  '+c.repaint[j].oldColor+' '+c.repaint[j].element+' = '+ c.repaint[j].newColor;
      if lowercase(s) <> lowercase(s1) then differs:=true;
    end;
  end;
  //if yes, write it
  if differs then
  begin
    writeln(F,'repaint     = '+IntToStr(length(Processes[i].repaint)));
    for j:=0 to length(Processes[i].repaint)-1 do
    begin
      writeln(F,'  '+Processes[i].repaint[j].oldColor+' '+Processes[i].repaint[j].element+' = '+ Processes[i].repaint[j].newColor);
    end;
  end;
  if Processes[i].feedingEfficiency <> c.feedingEfficiency then
    writeln(F,'feedingEfficiency = '+Processes[i].feedingEfficiency);
  if Processes[i].vertLoc <> c.vertLoc then
    writeln(F,'vertLoc     = '+IntToStrVertLoc(Processes[i].vertLoc));
  if Processes[i].isOutput <> c.isOutput then
    writeln(F,'isOutput    = '+IntToStr(Processes[i].isOutput));
  if Processes[i].isActive <> c.isActive then
    writeln(F,'isActive    = '+IntToStr(Processes[i].isActive));
  if Processes[i].processType <> c.processType then
    writeln(F,'processType = '+Processes[i].processType);
  //check if limitation block needs to be written
  if length(Processes[i].limitations) <> length(c.limitations) then differs:=true
  else
  begin
    differs:=false;
    for j:=0 to length(Processes[i].limitations)-1 do
    begin
      if processes[i].limitations[j].tracerIsSmall=0 then
      begin
        if processes[i].limitations[j].elseProcess <> '' then
          s:='limitation  = '
                    +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                    +limitations[processes[i].limitations[j].limitationNum].tracer+' > '
                    +limitations[processes[i].limitations[j].limitationNum].value+' else '
                    +processes[i].limitations[j].elseProcess
        else
          s:='limitation  = '
                    +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                    +limitations[processes[i].limitations[j].limitationNum].tracer+' > '
                    +limitations[processes[i].limitations[j].limitationNum].value;
      end
      else  //tracerIsSmall=1
      begin
        if processes[i].limitations[j].elseProcess <> '' then
          s:='limitation  = '
                    +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                    +limitations[processes[i].limitations[j].limitationNum].tracer+' < '
                    +limitations[processes[i].limitations[j].limitationNum].value+' else '
                    +processes[i].limitations[j].elseProcess
        else
          s:='limitation  = '
                    +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                    +limitations[processes[i].limitations[j].limitationNum].tracer+' < '
                    +limitations[processes[i].limitations[j].limitationNum].value;
      end;
      if c.limitations[j].tracerIsSmall=0 then
      begin
        if c.limitations[j].elseProcess <> '' then
          s1:='limitation  = '
                    +ref_limitations[c.limitations[j].limitationNum].limitationType+' '
                    +ref_limitations[c.limitations[j].limitationNum].tracer+' > '
                    +ref_limitations[c.limitations[j].limitationNum].value+' else '
                    +c.limitations[j].elseProcess
        else
          s1:='limitation  = '
                    +ref_limitations[c.limitations[j].limitationNum].limitationType+' '
                    +ref_limitations[c.limitations[j].limitationNum].tracer+' > '
                    +ref_limitations[c.limitations[j].limitationNum].value;
      end
      else  //tracerIsSmall=1
      begin
        if c.limitations[j].elseProcess <> '' then
          s1:='limitation  = '
                    +ref_limitations[c.limitations[j].limitationNum].limitationType+' '
                    +ref_limitations[c.limitations[j].limitationNum].tracer+' < '
                    +ref_limitations[c.limitations[j].limitationNum].value+' else '
                    +c.limitations[j].elseProcess
        else
          s1:='limitation  = '
                    +ref_limitations[c.limitations[j].limitationNum].limitationType+' '
                    +ref_limitations[c.limitations[j].limitationNum].tracer+' < '
                    +ref_limitations[c.limitations[j].limitationNum].value;
      end;
      if lowercase(s) <> lowercase(s1) then differs:=true;
    end;
  end;
  //if yes, write it
  if differs then
  begin
    for j:=0 to length(processes[i].limitations)-1 do

    begin
      if processes[i].limitations[j].tracerIsSmall=0 then
      begin
        if processes[i].limitations[j].elseProcess <> '' then
          writeln(F,'limitation  = '
                    +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                    +limitations[processes[i].limitations[j].limitationNum].tracer+' > '
                    +limitations[processes[i].limitations[j].limitationNum].value+' else '
                    +processes[i].limitations[j].elseProcess)
        else
          writeln(F,'limitation  = '
                    +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                    +limitations[processes[i].limitations[j].limitationNum].tracer+' > '
                    +limitations[processes[i].limitations[j].limitationNum].value);
      end
      else  //tracerIsSmall=1
      begin
        if processes[i].limitations[j].elseProcess <> '' then
          writeln(F,'limitation  = '
                    +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                    +limitations[processes[i].limitations[j].limitationNum].tracer+' < '
                    +limitations[processes[i].limitations[j].limitationNum].value+' else '
                    +processes[i].limitations[j].elseProcess)
        else
          writeln(F,'limitation  = '
                    +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                    +limitations[processes[i].limitations[j].limitationNum].tracer+' < '
                    +limitations[processes[i].limitations[j].limitationNum].value);
      end;
    end;
    if length(processes[i].limitations)=0 then writeln(F,'limitation = none');
  end;
end;

procedure SaveSingleProcess(var F: TextFile; i: Integer);
var c: TErgProcess;
    s: String;
    j: Integer;
begin
      initErgProcess(c);
      writeln(F,'name        = '+Processes[i].name);
      writeln(F,'description = '+Processes[i].description);
      writeln(F,'turnover    = '+Processes[i].turnover);
      if Processes[i].comment <> c.comment then
        writeln(F,'comment     = '+Processes[i].comment);

      //now, create equation string
      s:=GetInputOutputEquation(Processes[i]);
      writeln(F,'equation    = '+s);
      if length(processes[i].repaint)>0 then
      begin
        writeln(F,'repaint     = '+IntToStr(length(Processes[i].repaint)));
        for j:=0 to length(Processes[i].repaint)-1 do
        begin
          writeln(F,'  '+Processes[i].repaint[j].oldColor+' '+Processes[i].repaint[j].element+' = '+ Processes[i].repaint[j].newColor);
        end;
      end;
      if Processes[i].feedingEfficiency <> c.feedingEfficiency then
        writeln(F,'feedingEfficiency = '+Processes[i].feedingEfficiency);
      if Processes[i].vertLoc <> c.vertLoc then
        writeln(F,'vertLoc     = '+IntToStrVertLoc(Processes[i].vertLoc));
      if Processes[i].isOutput <> c.isOutput then
        writeln(F,'isOutput    = '+IntToStr(Processes[i].isOutput));
      if Processes[i].isActive <> c.isActive then
        writeln(F,'isActive    = '+IntToStr(Processes[i].isActive));
      if Processes[i].processType <> c.processType then
        writeln(F,'processType = '+Processes[i].processType);
      for j:=0 to length(processes[i].limitations)-1 do
      begin
        if processes[i].limitations[j].tracerIsSmall=0 then
        begin
          if processes[i].limitations[j].elseProcess <> '' then
            writeln(F,'limitation  = '
                      +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                      +limitations[processes[i].limitations[j].limitationNum].tracer+' > '
                      +limitations[processes[i].limitations[j].limitationNum].value+' else '
                      +processes[i].limitations[j].elseProcess)
          else
            writeln(F,'limitation  = '
                      +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                      +limitations[processes[i].limitations[j].limitationNum].tracer+' > '
                      +limitations[processes[i].limitations[j].limitationNum].value);
        end
        else  //tracerIsSmall=1
        begin
          if processes[i].limitations[j].elseProcess <> '' then
            writeln(F,'limitation  = '
                      +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                      +limitations[processes[i].limitations[j].limitationNum].tracer+' < '
                      +limitations[processes[i].limitations[j].limitationNum].value+' else '
                      +processes[i].limitations[j].elseProcess)
          else
            writeln(F,'limitation  = '
                      +limitations[processes[i].limitations[j].limitationNum].limitationType+' '
                      +limitations[processes[i].limitations[j].limitationNum].tracer+' < '
                      +limitations[processes[i].limitations[j].limitationNum].value);
        end;
      end;

  //Save extra tags
  SaveExtraTagsRecord(processes[i].ExtraTagsRecord, F, 12);
end;

procedure SaveProcesses(filename: String);
var F: TextFile;
    i: Integer;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  AssignFile(F,FileName);
  rewrite(F);
    WriteProcessesHeader(F);
    for i:=0 to length(Processes)-1 do
    begin
      SaveSingleProcess(F,i);
      writeln(F,'***********************');
    end;
  closefile(F);
end;

//************************** CElements **********************************//

function LoadSingleCElement(var F: TextFile): Integer; //returns 0 if succeeded
var
  s,s1,s2: String;
  i, p: Integer;
begin
  result:=0;
  setLength(CElements,length(CElements)+1);
  i:=length(CElements)-1;
  initErgCElement(CElements[length(CElements)-1]);
  s1:=''; s2:='';
  while not EOF(F) do
  begin
      readln(F,s);
      splitEquation(s,s1,s2);
      if s1 <> '' then
      begin
        if (copy(s1,1,1)='*') then            //CElement ended
          break
        else if lowercase(s1)='element' then
        begin
          CElements[i].element:=s2;
          for p:=0 to i-1 do                  //if celement already exists, update it
            if (lowercase(celements[i].element)=lowercase(celements[p].element)) and
              (lowercase(celements[i].color)=lowercase(celements[p].color)) then
            begin
              i:=p;
              setLength(celements, length(celements)-1);
            end;
        end
        else if lowercase(s1)='color' then
        begin
          CElements[length(CElements)-1].color:=s2;
          for p:=0 to i-1 do                  //if celement already exists, update it
            if (lowercase(celements[i].element)=lowercase(celements[p].element)) and
              (lowercase(celements[i].color)=lowercase(celements[p].color)) then
            begin
              i:=p;
              setLength(celements, length(celements)-1);
            end;
        end
        else if lowercase(s1)='description' then
          CElements[i].description:=s2
        else if lowercase(s1)='isaging' then
          CElements[i].isAging:=s2
        else if lowercase(s1)='atmosdep' then
          CElements[i].atmosdep:=StrToInt(s2)
        else if lowercase(s1)='riverdep' then
          CElements[i].riverdep:=StrToInt(s2)
        else if lowercase(s1)='istracer' then
          CElements[i].isTracer:=s2
        else if lowercase(s1)='longname_prefix' then
          CElements[i].longname_prefix:=s2
        else if lowercase(s1)='longname_suffix' then
          CElements[i].longname_suffix:=s2
        else if lowercase(s1)='duplicatedtracers' then
          CElements[i].duplicatedTracers:=s2
        else if lowercase(s1)='comment' then
          CElements[i].comment:=s2
        else //undefined variables were set, add as extra tag
          AddExtraTagFromText(s1, s2, CElements[i].ExtraTagsRecord);
        //annoying fix for that elements also must have names
        CElements[i].name := CElements[i].color+'_'+CElements[i].element;
      end;
  end; //while ended: now check if a CElement with this name exists
  if CElements[length(CElements)-1].element = '' then
    setLength(CElements,length(CElements)-1);
end;

function LoadCElements(filename: String; loadAddOn: Boolean=false):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  if not loadAddOn then setLength(CElements,0);
  AssignFile(F,filename);
  reset(F);
  result:=0;
    while not EOF(F) do
      result:=max(LoadSingleCElement(F),result);
  closefile(F);
end;

procedure WriteCElementsHeader(var F: TextFile);
begin
  writeln(F,'! BioCGT CElements file');
  writeln(F,'! *********************');
  writeln(F,'! properties of CElements: (colored elements to be traced)');
  writeln(F,'! element=           internal name of element, e.g., "N"');
  writeln(F,'! color=             e.g. "red", may not contain spaces');
  writeln(F,'! description=       e.g. "nitrogen from Oder river", default=""');
  writeln(F,'! atmosDep=          1=atmospheric deposition of marked tracers may occur, 0=not (default)');
  writeln(F,'! isAging=           1=accumulates time since entering the system, 0=does not (default)');
  writeln(F,'! isTracer=          1=store total Element content in a separate tracer, 0=do not (default)');
  writeln(F,'!                    setting isAging=1 implies isTracer=1');
  writeln(F,'! riverDep=          1=river deposition of marked tracers may occur, 0=not (default)');
  writeln(F,'! longname_prefix=   prefix to be prepended to a tracers longname when painted; only applied');
  writeln(F,'!                    when tracer-based naming convention is acitve (default); default: ""');
  writeln(F,'! longname_suffix    suffix to be appended to a tracers longname when painted; only applied');
  writeln(F,'!                    when tracer-based naming convention is acitve (default);');
  writeln(F,'!                    default: "" (longname gets same suffix as name)');
  writeln(F,'! duplicatedTracers= "all" (default) or a semicolon-separated list of tracer names that should be duplicated');
  writeln(F,'! comment=           any comments');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleCElementRef(var F: TextFile; i: Integer; c: TErgCElement; AlwaysShowDefault: Boolean=true);
begin
  writeln(F,'element           = '+CElements[i].element);
  writeln(F,'color             = '+CElements[i].color);
  if (CElements[i].description <> c.description) or AlwaysShowDefault then
    writeln(F,'description       = '+CElements[i].description);
  if CElements[i].isAging <> c.isAging then
    writeln(F,'isAging           = '+CElements[i].isAging);
  if CElements[i].atmosDep <> c.atmosDep then
    writeln(F,'atmosDep          = '+IntToStr(CElements[i].atmosDep));
  if CElements[i].riverDep <> c.riverDep then
    writeln(F,'riverDep          = '+IntToStr(CElements[i].riverDep));
  if CElements[i].isTracer <> c.isTracer then
    writeln(F,'isTracer          = '+CElements[i].isTracer);
  if CElements[i].longname_prefix <> c.longname_prefix then
    writeln(F,'longname_prefix   = '+CElements[i].longname_prefix);
  if CElements[i].longname_suffix <> c.longname_suffix then
    writeln(F,'longname_suffix   = '+CElements[i].longname_suffix);
  if CElements[i].duplicatedTracers <> c.duplicatedTracers then
    writeln(F,'duplicatedTracers = '+CElements[i].duplicatedTracers);
  if CElements[i].comment <> c.comment then
    writeln(F,'comment           = '+CElements[i].comment);

  //Save extra tags
  SaveExtraTagsRecord(CElements[i].ExtraTagsRecord, F, 18);
end;

procedure SaveSingleCElement(var F: TextFile; i: Integer);
var c: TErgCElement;
begin
  InitErgCElement(c);
  SaveSingleCElementRef(F,i,c);
end;

procedure SaveCElements(filename: String);
var F: TextFile;
    i: Integer;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  AssignFile(F,FileName);
  rewrite(F);
    WriteCElementsHeader(F);
    for i:=0 to length(CElements)-1 do
    begin
      SaveSingleCElement(F,i);
      writeln(F,'***********************');
    end;
  closefile(F);
end;

//************************** modelInfos **********************************//

function LoadModelInfos(filename: String):Integer; //returns 0 if succeeded
var F: TextFile;
    s,s1,s2: String;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  result:=1;
  initErgModelInfo(ModelInfos);
  AssignFile(F,filename);
  reset(F);
    s1:=''; s2:='';
    while not EOF(F) do
    begin
      readln(F,s);
      while copy(trim(s),1,1)='*' do readln(F,s);
      splitEquation(s,s1,s2);
      if s1 <> '' then
      begin
        if lowercase(s1)='name' then
          ModelInfos.name:=s2
        else if lowercase(s1)='version' then
          ModelInfos.version:=s2
        else if lowercase(s1)='description' then
          ModelInfos.description:=s2
        else if lowercase(s1)='author' then
          ModelInfos.author:=s2
        else if lowercase(s1)='contact' then
          ModelInfos.contact:=s2
        else if lowercase(s1)='templatepath' then
          if (copy(s2,length(s2),1)='\') or
             (copy(s2,length(s2),1)='/') then
            ModelInfos.templatePath:=s2
          else
            ModelInfos.templatePath:=s2+'/'
        else if lowercase(s1)='outputpath' then
          if (copy(s2,length(s2),1)='\') or
             (copy(s2,length(s2),1)='/') then
            ModelInfos.outputPath:=s2
          else
            ModelInfos.outputPath:=s2+'/'
        else if lowercase(s1)='ageepsilon' then
          ModelInfos.ageEpsilon:=s2
        else if lowercase(s1)='autolimitprocesses' then
          ModelInfos.autoLimitProcesses:=StrToInt(s2)
        else if lowercase(s1)='automassclassprop' then
          ModelInfos.autoMassClassProp:=StrToInt(s2)
        else if lowercase(s1)='autosplitcolors' then
          ModelInfos.autoSplitColors:=StrToInt(s2)
        else if lowercase(s1)='autosortmoving' then
          ModelInfos.autoSortMoving:=StrToInt(s2)
        else if lowercase(s1)='debugmode' then
          Modelinfos.debugMode:=StrToInt(s2)
        else if lowercase(s1)='autowrapf' then
          ModelInfos.autoWrapF:=StrToInt(s2)
        else if lowercase(s1)='autowrapf90' then
          ModelInfos.autoWrapF90:=StrToInt(s2)
        else if lowercase(s1)='realsuffixf90' then
          ModelInfos.RealSuffixF90:=s2
        else if lowercase(s1)='autounixoutput' then
          ModelInfos.autoUnixOutput:=StrToInt(s2)
        else if lowercase(s1)='inactiveprocesstypes' then
          ModelInfos.inactiveProcessTypes:=s2
        else if lowercase(s1)='numphysicalparametersvarying' then
          ModelInfos.numPhysicalParametersVarying:=StrToInt(s2)
        else if lowercase(s1)='conversionfactorsed' then
          ModelInfos.conversionFactorSED:=s2
        else if lowercase(s1)='conversionfactorsur' then
          ModelInfos.conversionFactorSUR:=s2
        else if lowercase(s1)='comment' then
          ModelInfos.comment:=s2
        else //undefined variables were set, add as extra tag
          AddExtraTagFromText(s1, s2, ModelInfos.ExtraTagsRecord);
      end;
    end;
  closefile(F);
  if result=1 then result:=0;
end;

procedure WriteModelInfosHeader(var F: TextFile);
begin
  writeln(F,'! BioCGT ModelInfos file');
  writeln(F,'! **********************');
  writeln(F,'! properties of ModelInfos:');
  writeln(F,'! name=                 bio-model short name or abbreviation');
  writeln(F,'! description=          bio-model long name');
  writeln(F,'! version=              bio-model version');
  writeln(F,'! author=               bio-model author(s)');
  writeln(F,'! contact=              e.g. e-mail adress of bio-model author');
  writeln(F,'! ageEpsilon=           small value used for preventing zero division for age calculation; default="1.0e-20"');
  writeln(F,'! autoBurialFluxes=     1=auto-generate fluxes for burial of colored elements with isTracer=1; 0=do not (default)');
  writeln(F,'! autoLimitProcesses=   1=add limitations to all processes that stop them when one of their precursors with isPositive=1 becomes zero (default); 0=do not');
  writeln(F,'! autoMassClassProp=    0=manual mass-class propagation processes (default); 1=mass-class propagation when upper mass limit is reached; 2=advanced propagation; 3=age-class propagation at beginning of each year');
  writeln(F,'! autoSortMoving=       1=sort tracers: not vertically moving first, then vertically moving; 0=do not (default)');
  writeln(F,'! autoSplitColors=      1=split tracers and processes according to colored elements (default); 0=do not');
  writeln(F,'! autoUnixOutput=       1=enforce Unix line-feed output on Windows systems; 0=do not (default)');
  writeln(F,'! autoWrapF=            1=auto-wrap too long lines in all files with ".f" or ".F" extension (default); 0=do not');
  writeln(F,'! autoWrapF90=          1=auto-wrap too long lines in all files with ".f90" or ".F90" extension; 0=do not (default)');
  writeln(F,'! conversionFactorSED=  factor to convert between vertLoc=WAT and vertLoc=SED tracers, default="(cgt_cellheight*cgt_density)"');
  writeln(F,'! conversionFactorSUR=  factor to convert between vertLoc=WAT and vertLoc=SUR tracers, default="(cgt_cellheight*cgt_density)"');
  writeln(F,'! realSuffixF90=        e.g. "8", appends a suffix (e.g. _8) to all real values in .f90 files which do not yet contain a suffix (do not include the underscore here); default='' meaning no suffix is added.');
  writeln(F,'! debugMode=            1=debug mode (output of all values); 0=output only of those values with output=1 (default)');
  writeln(F,'! inactiveProcessTypes= semicolon-separated list of process types that are set inactive, e.g. because they are represented in the host model, e.g. "gas_exchange; sedimentation; resuspension"');
  writeln(F,'! numPhysicalParametersVarying= number of physical parameters that are allowed to vary in an inverse-modelling approach (default=0)');
  writeln(F,'! outputPath=           path where to write the output files');
  writeln(F,'! templatePath=         path to the code template files');
  writeln(F,'! comment=              comments about the current model version');
  writeln(F,'!');
  writeln(F,'! use ! for comments');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveModelInfos(filename: String);
var F: TextFile;
    c: TErgModelInfo;
begin
  DefaultFormatSettings.DecimalSeparator:='.';
  InitErgModelInfo(c);
  AssignFile(F,FileName);
  rewrite(F);
  WriteModelInfosHeader(F);
  writeln(F,'name              = '+ModelInfos.name);
  writeln(F,'version           = '+ModelInfos.version);
  writeln(F,'description       = '+ModelInfos.description);
  writeln(F,'author            = '+ModelInfos.author);
  writeln(F,'contact           = '+ModelInfos.contact);
  writeln(F,'templatePath      = '+ModelInfos.templatePath);
  writeln(F,'outputPath        = '+ModelInfos.outputPath);
  if ModelInfos.ageEpsilon <> c.ageEpsilon then
    writeln(F,'ageEpsilon        = '+ModelInfos.ageEpsilon);
  if ModelInfos.autoLimitProcesses<> c.autoLimitProcesses then
    writeln(F,'autoLimitProcesses = '+IntToStr(ModelInfos.autoLimitProcesses));
  if ModelInfos.autoMassClassProp <> c.autoMassClassProp then
    writeln(F,'autoMassClassProp = '+IntToStr(ModelInfos.autoMassClassProp));
  if ModelInfos.autoSortMoving <> c.autoSortMoving then
    writeln(F,'autoSortMoving    = '+IntToStr(ModelInfos.autoSortMoving));
  if ModelInfos.autoSplitColors <> c.autoSplitColors then
    writeln(F,'autoSplitColors   = '+IntToStr(ModelInfos.autoSplitColors));
  if ModelInfos.debugMode <> c.debugMode then
    writeln(F,'debugMode         = '+IntToStr(ModelInfos.debugMode));
  if ModelInfos.autoWrapF <> c.autoWrapF then
    writeln(F,'autoWrapF         = '+IntToStr(ModelInfos.autoWrapF));
  if ModelInfos.autoWrapF90 <> c.autoWrapF90 then
    writeln(F,'autoWrapF90       = '+IntToStr(ModelInfos.autoWrapF90));
  if ModelInfos.RealSuffixF90 <> c.RealSuffixF90 then
    writeln(F,'realSuffixF90     = '+ModelInfos.RealSuffixF90);
  if ModelInfos.autoUnixOutput <> c.autoUnixOutput then
    writeln(F,'autoUnixOutput    = '+IntToStr(ModelInfos.autoUnixOutput));
  if ModelInfos.inactiveProcessTypes <> c.inactiveProcessTypes then
    writeln(F,'inactiveProcessTypes = '+ModelInfos.InactiveProcessTypes);
  if ModelInfos.numPhysicalParametersVarying <> c.numPhysicalParametersVarying then
    writeln(F,'numPhysicalParametersVarying = '+IntToStr(ModelInfos.numPhysicalParametersVarying));
  if ModelInfos.conversionFactorSED <> c.conversionFactorSED then
    writeln(F,'conversionFactorSED = '+ModelInfos.conversionFactorSED);
  if ModelInfos.conversionFactorSUR <> c.conversionFactorSUR then
    writeln(F,'conversionFactorSUR = '+ModelInfos.conversionFactorSUR);
  if ModelInfos.comment <> c.comment then
    writeln(F,'comment              = '+ModelInfos.comment);
  //Save extra tags
  SaveExtraTagsRecord(ModelInfos.ExtraTagsRecord, F, 18);
  closefile(F);
end;

//******** Index Generation **********************************************//

procedure SortMovingTracers;
//sort tracers: not vertically moving tracers first, then moving ones  (for MOM4 compatibility)
var
  tracer: tErgTracer;
  i, j: Integer;
begin
  tracer:=tErgTracer.create;
  for i:=0 to length(tracers)-1 do
  begin
    if (tracers[i].vertSpeed<>'0') and (tracers[i].vertSpeed<>'0.0') then //this is a moving tracer
      for j:=i+1 to length(tracers)-1 do
        if (tracers[j].vertSpeed='0') or (tracers[j].vertSpeed='0.0') then //this tracer does not move
        begin
          //exchange these two tracers
          copyErgTracer(tracers[i],tracer);
          copyErgTracer(tracers[j],tracers[i]);
          copyErgTracer(tracer,tracers[j]);
          break;
        end;
  end;
end;

function GenerateIndexes(sort_moving_notmoving: Boolean=false): String;
var i,j,k,l, ii: Integer;
    s: String;
    e: TParseExpression;
    precedingStuff: String;
    numerator, denominator: String;
    ratestring: String;
begin
  s:='';
  DefaultFormatSettings.DecimalSeparator:='.';

  if sort_moving_notmoving then SortMovingTracers;  //sort tracers: not vertically moving tracers first, then moving ones

  ConstantsList.Free; AuxiliariesList.Free; ProcessesList.Free;
  ConstantsList:=TStringList.Create;
  AuxiliariesList:=TStringList.Create;
  ProcessesList:=TStringList.Create;

  for i:=0 to length(constants)-1 do
  begin
    if pos('$', constants[i].name)<=0 then
    begin
      if pos(';',constants[i].valueString)<=0 then
        ConstantsList.Add(constants[i].name+'='+Constants[i].valueString);
      if constants[i].valueString='' then
        s:=s+'Constant '+constants[i].name+' has no value. If you use vector tracers, check if you forgot a semicolon!'+chr(13)
      else
      begin //use parser to get value
        e:=TParseExpression.create;
        e.text:=constants[i].name;
        e.Evaluate(ConstantsList);
        if e.error then s:=s+'Constant '+constants[i].name+': Error calculating the value: '+e.errorText+chr(13)
        else constants[i].value:=e.value;
        e.Free;
      end;
    end;
  end;

  for i:=0 to length(auxiliaries)-1 do
    if pos(';',auxiliaries[i].formula)<=0 then
      AuxiliariesList.Add(auxiliaries[i].name+'='+auxiliaries[i].formula);

  for i:=0 to length(processes)-1 do
    if pos(';',processes[i].turnover)<=0 then
      ProcessesList.Add(processes[i].name+'='+Processes[i].turnover);

  for i:=0 to length(tracers)-1 do
  begin
    setLength(tracers[i].myRateList,0);
    tracers[i].myChildOf:=-1;
    tracers[i].myTracerAbove:=-1;
    tracers[i].myTracerBelow:=-1;
    for j:=0 to length(tracers)-1 do
    begin
      if lowercase(tracers[j].name) = lowercase(tracers[i].childOf) then
        tracers[i].myChildOf:=j;
      if lowercase(tracers[j].name) = lowercase(tracers[i].tracerAbove) then
        tracers[i].myTracerAbove:=j;
      if lowercase(tracers[j].name) = lowercase(tracers[i].tracerBelow) then
        tracers[i].myTracerBelow:=j;
    end;
    if (tracers[i].myChildOf=-1) and (lowercase(tracers[i].childOf)<>'none') then
      s:=s+'Tracer "'+tracers[i].name+'" is child of "'+tracers[i].childOf+'", but that tracer does not exist.'+chr(13);
    if (tracers[i].myTracerAbove=-1) and (lowercase(tracers[i].tracerAbove)<>'none') then
      s:=s+'Tracer "'+tracers[i].name+'" is below "'+tracers[i].tracerAbove+'" (=tracerAbove), but that tracer does not exist.'+chr(13);
    if (tracers[i].myTracerBelow=-1) and (lowercase(tracers[i].tracerBelow)<>'none') then
      s:=s+'Tracer "'+tracers[i].name+'" is above "'+tracers[i].tracerBelow+'" (=tracerBelow), but that tracer does not exist.'+chr(13);
    if pos(copy(tracers[i].vertSpeed,1,1),'-0123456789.')>0 then
    begin
      try
        tracers[i].myVertspeedValue:=StrToFloat(tracers[i].vertSpeed);
      except
        on eConvertError do
        begin
          tracers[i].myVertspeedValue:=0;
          for j:=0 to length(constants)-1 do
            if lowercase(tracers[i].vertSpeed)=lowercase(constants[j].name) then
              tracers[i].myVertspeedValue:=constants[j].value;
        end;
      end;
    end
    else
    begin
        tracers[i].myVertspeedValue:=0;
        for j:=0 to length(constants)-1 do
          if lowercase(tracers[i].vertSpeed)=lowercase(constants[j].name) then
            tracers[i].myVertspeedValue:=constants[j].value;
    end;
    if pos(copy(tracers[i].opacity,1,1),'-0123456789.')>0 then
    begin
      try
        tracers[i].myOpacityValue:=StrToFloat(tracers[i].opacity);
      except
        on eConvertError do
        begin
          tracers[i].myOpacityValue:=0;
          for j:=0 to length(constants)-1 do
            if lowercase(tracers[i].opacity)=lowercase(constants[j].name) then
              tracers[i].myOpacityValue:=constants[j].value;
        end;
      end;
    end
    else
    begin
        tracers[i].myOpacityValue:=0;
        for j:=0 to length(constants)-1 do
          if lowercase(tracers[i].opacity)=lowercase(constants[j].name) then
            tracers[i].myOpacityValue:=constants[j].value;
    end;
    for j:=0 to length(tracers[i].contents)-1 do
    begin
      tracers[i].contents[j].myElementNum:=-1;
      if tracers[i].isCombined=0 then  // the default case, tracer contains some elements
      begin
        for k:=0 to length(elements)-1 do
          if elements[k].name = tracers[i].contents[j].element then
            tracers[i].contents[j].myElementNum := k;
        if tracers[i].contents[j].myElementNum = -1 then
          s:=s+'Tracer "'+tracers[i].name+'" contains element "'+tracers[i].contents[j].element+'", but that element does not exist.'+chr(13);
      end
      else    // tracer is a combined tracer and contains no elements, but other virtual (isActive=false) tracers instead
      begin
        for k:=0 to length(tracers)-1 do
          if tracers[k].name = tracers[i].contents[j].element then
          begin
            tracers[i].contents[j].myElementNum := k;
            if tracers[k].isCombined>0 then
              s:=s+'Combined tracer "'+tracers[i].name+'" contains tracer "'+tracers[i].contents[j].element+'", but that is also combined and not virtual.'+chr(13);
          end;
        if tracers[i].contents[j].myElementNum = -1 then
          s:=s+'Combined tracer "'+tracers[i].name+'" contains tracer "'+tracers[i].contents[j].element+'", but that tracer does not exist.'+chr(13);
      end;

      e:=TParseExpression.Create;
      e.text:=tracers[i].contents[j].amount;
      e.Evaluate(ConstantsList);
      if e.error then
      begin
        s:=s+'Tracer "'+tracers[i].name+'" contains element "'+tracers[i].contents[j].element+'" in amount "'+tracers[i].contents[j].amount+'", which cannot be evaluated, error: '+e.errorText+chr(13);
        tracers[i].contents[j].myAmount:=0;
      end
      else
      begin
        tracers[i].contents[j].myAmount:=e.value;
      end;
      e.Free;
    end;
  end;

  for i:=0 to length(limitations)-1 do
  begin
    limitations[i].myTracerNum:=-1;
    for j:=0 to length(tracers)-1 do
      if lowercase(limitations[i].tracer)=lowercase(tracers[j].name) then
        limitations[i].myTracerNum:=j;
    if limitations[i].myTracerNum=-1 then
      s:=s+'A process limitation exists with tracer "'+limitations[i].tracer+', which does not exist.'+chr(13);
  end;

  for i:=0 to length(processes)-1 do
  begin
    processes[i].myIsStiff:=0;
    processes[i].myStiffTracerNum:=-1;
    for j:=0 to length(processes[i].input)-1 do
    begin
      //find numbers of tracers which are consumed
      processes[i].input[j].myTracerNum:=-1;
      for k:=0 to length(tracers)-1 do
      begin
        if lowercase(processes[i].input[j].tracer) = lowercase(tracers[k].name) then
          processes[i].input[j].myTracerNum:=k;
        if tracers[k].dimension>0 then
        begin
          if copy(lowercase(processes[i].input[j].tracer),1,length(tracers[k].name)+2)=lowercase(tracers[k].name)+'_$' then  // check if output of a vector tracer occurs
            processes[i].input[j].myTracerNum:=k
          else
            for l:=1 to tracers[k].dimension do
            begin
              if (lowercase(processes[i].input[j].tracer) = lowercase(tracers[k].name)+'_'+IntToStr(l)) and (tracers[k].massLimits='') then processes[i].input[j].myTracerNum:=k;
              if (lowercase(processes[i].input[j].tracer) = lowercase(tracers[k].name)+'_'+IntToStr(l)+'_mass') and (tracers[k].massLimits<>'') then processes[i].input[j].myTracerNum:=k;
              if (lowercase(processes[i].input[j].tracer) = lowercase(tracers[k].name)+'_'+IntToStr(l)+'_abd') and (tracers[k].massLimits<>'') then processes[i].input[j].myTracerNum:=k;
            end;
        end;
      end;
      if processes[i].input[j].myTracerNum=-1 then
        if (pos('$',processes[i].input[j].tracer)<=0) then //prevent error message being displayed if input is a vector tracer
          s:=s+'Process "'+processes[i].name+'" has input of tracer "'+processes[i].input[j].tracer+'", which is not defined.'+chr(13);

      //check if the amount of all consumed tracers is properly defined.
      e:=TParseExpression.Create;
      e.text:=processes[i].input[j].amount;
      e.Evaluate(ConstantsList);
      if e.error then
      begin
        s:=s+'Process "'+processes[i].name+'" has input of tracer "'+processes[i].input[j].tracer+'" in amount "'+processes[i].input[j].amount+'", which cannot be evaluated, error: '+e.errorText+chr(13);
        processes[i].input[j].myAmount:=0;
      end
      else
      begin
        processes[i].input[j].myAmount:=e.value;
      end;
      e.Free;

      //check if stiff tracer is consumed
      k:=processes[i].input[j].myTracerNum;
      if k<>-1 then
      begin
        if tracers[k].isStiff<>0 then
        begin
          if processes[i].myIsStiff<>0 then
            s:=s+'Process "'+processes[i].name+'" has input of more than one tracer with isStiff/=0.'+chr(13)
          else
          begin
            processes[i].myIsStiff:=tracers[k].isStiff;
            processes[i].myStiffTracerNum:=k;
          end;
        end;
      end;
    end;
    for j:=0 to length(processes[i].output)-1 do
    begin
      processes[i].output[j].myTracerNum:=-1;
      for k:=0 to length(tracers)-1 do
      begin
        if lowercase(processes[i].output[j].tracer) = lowercase(tracers[k].name) then
          processes[i].output[j].myTracerNum:=k;
        if tracers[k].dimension>0 then
        begin
          if copy(lowercase(processes[i].output[j].tracer),1,length(tracers[k].name)+2)=lowercase(tracers[k].name)+'_$' then  // check if output of a vector tracer occurs
            processes[i].output[j].myTracerNum:=k
          else // maybe only one instance of the vector tracer occurs in output, check that
            for l:=1 to tracers[k].dimension do
            begin
              if (lowercase(processes[i].output[j].tracer) = lowercase(tracers[k].name)+'_'+IntToStr(l)) and (tracers[k].massLimits='') then processes[i].output[j].myTracerNum:=k;
              if (lowercase(processes[i].output[j].tracer) = lowercase(tracers[k].name)+'_'+IntToStr(l)+'_mass') and (tracers[k].massLimits<>'') then processes[i].output[j].myTracerNum:=k;
              if (lowercase(processes[i].output[j].tracer) = lowercase(tracers[k].name)+'_'+IntToStr(l)+'_abd') and (tracers[k].massLimits<>'') then processes[i].output[j].myTracerNum:=k;
            end;
        end;
      end;
      if processes[i].output[j].myTracerNum=-1 then
        if (pos('$',processes[i].output[j].tracer)<=0) then //prevent error message being displayed if output is a stage-resolving vector tracer
          s:=s+'Process "'+processes[i].name+'" has output of tracer "'+processes[i].output[j].tracer+'", which is not defined.'+chr(13);

      e:=TParseExpression.Create;
      e.text:=processes[i].output[j].amount;
      e.Evaluate(ConstantsList);
      if e.error then
      begin
        s:=s+'Process "'+processes[i].name+'" has output of tracer "'+processes[i].output[j].tracer+'" in amount "'+processes[i].output[j].amount+'", which cannot be evaluated, error: '+e.errorText+chr(13);
        processes[i].output[j].myAmount:=0;
      end
      else
      begin
        processes[i].output[j].myAmount:=e.value;
      end;
      e.Free;

      //check if stiff tracer is produced
      k:=processes[i].output[j].myTracerNum;
      if k<>-1 then
      begin
        if tracers[k].isStiff<>0 then
        begin
          if processes[i].myIsStiff<>0 then
            s:=s+'Process "'+processes[i].name+'" consumes and produces a tracer with isStiff/=0, or produces more than one tracer with isStiff=3.'+chr(13)
          else if tracers[k].isStiff=3 then
          begin
            processes[i].myIsStiff:=tracers[k].isStiff;
            processes[i].myStiffTracerNum:=k;
          end;
        end;
      end;
    end;
    for j:=0 to length(processes[i].repaint)-1 do
    begin
      processes[i].repaint[j].myElementNum:=-1;
      for k:=0 to length(elements)-1 do
        if lowercase(processes[i].repaint[j].element) = lowercase(elements[k].name) then
          processes[i].repaint[j].myElementNum:=k;
      if lowercase(processes[i].repaint[j].element)='all' then
        processes[i].repaint[j].myElementNum:=-2;
      if processes[i].repaint[j].myElementNum=-1 then
        s:=s+'Process "'+processes[i].name+'" repaints element "'+processes[i].repaint[j].element+'", which is not defined.'+chr(13);
    end;
    for j:=0 to length(processes[i].limitations)-1 do
    begin
      processes[i].limitations[j].myElseProcessNum:=-1;
      //elseProcess may be of the shape  "anything * processName", where anything provides a ratio (speed of replacement process)/(speed of original process)
      for k:=i+1 to length(processes)-1 do  //the process in "elseProcess" must occur behind this process to avoid cyclic definition
        if length(trim(processes[i].limitations[j].elseProcess))>=length(processes[k].name) then
          if copy(trim(lowercase(processes[i].limitations[j].elseProcess)),
                  length(trim(processes[i].limitations[j].elseProcess))-length(processes[k].name)+1,
                  length(processes[k].name)) = lowercase(processes[k].name) then //the process name is found at the end
          begin
            if length(trim(processes[i].limitations[j].elseProcess))=length(processes[k].name) then //elseProcess = processName
            begin
              processes[i].limitations[j].myElseProcessNum:=k;
              processes[i].limitations[j].myElseProcessRatio:='1.0';
            end
            else  //something precedes the processName, we hope it is a factor
            begin
              precedingStuff:=trim(copy(trim(processes[i].limitations[j].elseProcess),1,length(trim(processes[i].limitations[j].elseProcess))-length(processes[k].name)));
              if copy(precedingStuff,length(precedingStuff),1)='*' then //the preceding stuff is really a factor
              begin
                processes[i].limitations[j].myElseProcessNum:=k;
                processes[i].limitations[j].myElseProcessRatio:='('+copy(precedingStuff,1,length(precedingStuff)-1)+')';
              end;
            end;
          end;
      if (processes[i].limitations[j].myElseProcessNum=-1) and (processes[i].limitations[j].elseProcess <> '') then
        s:=s+'Process "'+processes[i].name+'" , if limited by tracer "'+limitations[processes[i].limitations[j].limitationNum].tracer+'", shall be replaced by process "'+processes[i].limitations[j].elseProcess+'" which is neither a process name which occurs behind "'+processes[i].name+'" in the processes list, nor such a process name preceded by a factor.'+chr(13);
    end;
  end;

  //another process loop for determining the patankar factors
  for i:=0 to length(processes)-1 do
    if processes[i].myIsStiff=0 then
      processes[i].myStiffFactor:='1.0'
    else
    begin
      k:=processes[i].myStiffTracerNum;
      numerator:=tracers[k].name+' + cgt_timestep*(0.0';
      for ii:=0 to length(processes)-1 do
      begin
        if processes[ii].isActive=1 then
        begin
          for j:=0 to length(processes[ii].output)-1 do   //all products of this process are checked
          begin
            if processes[ii].output[j].myTracerNum=k then   //the process i produces tracer t
            begin
              if (processes[ii].output[j].amount='1') or (processes[ii].output[j].amount='1.0') then
                ratestring:='+ '+processes[ii].name
              else
                ratestring:='+ ('+processes[ii].name+')*('+processes[ii].output[j].amount+')';
              //flat processes have rates in [mmol/m**2/day)], non-flat tracers have unit [mmol/m**3]
              if (tracers[k].vertLoc=0) and (processes[ii].vertLoc=1) then
                ratestring:=ratestring+'/'+modelinfos.conversionFactorSED;
              if (tracers[k].vertLoc=0) and (processes[ii].vertLoc=2) then
                ratestring:=ratestring+'/'+modelinfos.conversionFactorSUR;
              if (tracers[k].vertLoc=3) and (processes[ii].vertLoc=0) then
                ratestring:=ratestring+'*'+modelinfos.conversionFactorSUR;
              numerator:=numerator+ratestring;
            end;
          end;
        end;
      end;
      numerator:=numerator+')';

      denominator:=tracers[k].name+' + cgt_timestep*(0.0';
      for ii:=0 to length(processes)-1 do
      begin
        if processes[ii].isActive=1 then
        begin
          for j:=0 to length(processes[ii].input)-1 do   //all products of this process are checked
          begin
            if processes[ii].input[j].myTracerNum=k then   //the process i consumes tracer t
            begin
              if (processes[ii].input[j].amount='1') or (processes[ii].input[j].amount='1.0') then
                ratestring:='+ '+processes[ii].name
              else
                ratestring:='+ ('+processes[ii].name+')*('+processes[ii].input[j].amount+')';
              //flat processes have rates in [mmol/m**2/day)], non-flat tracers have unit [mmol/m**3]
              if (tracers[k].vertLoc=0) and (processes[ii].vertLoc=1) then
                ratestring:=ratestring+'/'+modelinfos.conversionFactorSED;
              if (tracers[k].vertLoc=0) and (processes[ii].vertLoc=2) then
                ratestring:=ratestring+'/'+modelinfos.conversionFactorSUR;
              if (tracers[k].vertLoc=3) and (processes[ii].vertLoc=0) then
                ratestring:=ratestring+'*'+modelinfos.conversionFactorSUR;
              denominator:=denominator+ratestring;
            end;
          end;
        end;
      end;
      denominator:=denominator+')';
      denominator:='max('+denominator+',1.0e-30)';
      processes[i].myStiffFactor:='('+numerator+')/('+denominator+')';
      if processes[i].myIsStiff=1 then processes[i].myStiffFactor:='min('+processes[i].myStiffFactor+',1.0)';
    end;

  for i:=0 to length(cElements)-1 do
  begin
    cElements[i].myIsAging:=StrToInt(cElements[i].isAging);
    if cElements[i].myIsAging=1 then cElements[i].isTracer:='1';
    cElements[i].myIsTracer:=StrToInt(cElements[i].isTracer);
    cElements[i].myElementNum:=-1;
    for j:=0 to length(elements)-1 do
      if lowercase(cElements[i].element) = lowercase(elements[j].name) then
        cElements[i].myElementNum:=j;
    if cElements[i].myElementNum=-1 then
      s:=s+'Colored elements contain element "'+cElements[i].element+'", which is not defined.'+chr(13);
  end;

  //now check which auxiliaries and tracers need to be calculated before the zIntegrals
  for j:=0 to length(auxiliaries)-1 do
    auxiliaries[j].myCalcBeforeZIntegral:=0;
  //First possibility: Auxiliary variable i is zIntegral over auxiliary variable j => j needs to be precalculated
  for i:=length(auxiliaries)-1 downto 0 do
    if ((auxiliaries[i].isZIntegral=1) and (auxiliaries[i].calcAfterProcesses=0)) then
    begin
      for j:=0 to i-1 do
      begin
        if    ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[1])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[3])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[4])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[5])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[6])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[7])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[8])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[9])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].formula) then
           auxiliaries[j].myCalcBeforeZIntegral:=1;
      end;
    end;
  //second possibility: Auxiliary variable i has vertLoc=WAT, but depends on j with vertLoc=SED or SUR => j needs to be precalculated
  //exception: j has isZIntegral, in this case, it is automatically precalculated
  for i:=length(auxiliaries)-1 downto 0 do
    if (auxiliaries[i].vertLoc=0) then
    begin
      for j:=0 to i-1 do
        if ((auxiliaries[j].vertLoc=1) or (auxiliaries[j].vertLoc=2)) and (auxiliaries[j].isZIntegral=0) then
        begin
          if    ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[1])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[3])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[4])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[5])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[6])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[7])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[8])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[9])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].formula) then
           auxiliaries[j].myCalcBeforeZIntegral:=1;
        end;
    end;
  //third possibility: Process i has vertLoc=WAT, but depends on auxiliary variable j with vertLoc=SED or SUR => j needs to be precalculated
  //exception: j has isZIntegral, in this case, it is automatically precalculated
  for i:=length(processes)-1 downto 0 do
    if (processes[i].vertLoc=0) then
    begin
      for j:=0 to length(auxiliaries)-1 do
        if ((auxiliaries[j].vertLoc=1) or (auxiliaries[j].vertLoc=2)) and (auxiliaries[j].isZIntegral=0) then
        begin
          if ContainedInString(auxiliaries[j].name,processes[i].turnover) then
            auxiliaries[j].myCalcBeforeZIntegral:=1;
        end;
    end;
  //Last possibility: Auxiliary variable i needs to be calculated before, but depends on j => j needs to be precalculated, too.
  //one exception is if i has vertLoc=SED and isZIntegral=0, but j has isZIntegral=1, then, j is automatically precalculated.
  for i:=length(auxiliaries)-1 downto 0 do
    if auxiliaries[i].myCalcBeforeZIntegral=1 then
    begin
      for j:=0 to i-1 do
      begin
        if    ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[1])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[3])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[4])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[5])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[6])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[7])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[8])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[9])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].formula) then
           if not ((auxiliaries[i].vertLoc=1) and (auxiliaries[i].isZIntegral=0) and (auxiliaries[j].isZIntegral=1)) then
             auxiliaries[j].myCalcBeforeZIntegral:=1;
      end;
    end;
  //Now we check which tracers we need to do the precalculation
  for j:=0 to length(tracers)-1 do
    tracers[j].myCalcBeforeZIntegral:=0;
  for i:=0 to length(auxiliaries)-1 do
    if ((auxiliaries[i].isZIntegral=1) and (auxiliaries[i].calcAfterProcesses=0)) or (auxiliaries[i].myCalcBeforeZIntegral=1) then
    begin
      for j:=0 to length(tracers)-1 do
      begin
        if    ContainedInString(tracers[j].name,auxiliaries[i].temp[1])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[2])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[3])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[4])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[5])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[6])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[7])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[8])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[9])
           or ContainedInString(tracers[j].name,auxiliaries[i].temp[2])
           or ContainedInString(tracers[j].name,auxiliaries[i].formula) then
           tracers[j].myCalcBeforeZIntegral:=1;
      end;
    end;
  //perform a consistency check: No variable with isZIntegral=1 may have myCalcBeforeZIntegral=1
  for i:=0 to length(auxiliaries)-1 do
    if (auxiliaries[i].isZIntegral=1) and (auxiliaries[i].myCalcBeforeZIntegral=1) then
      s:=s+'Auxiliary variable '+auxiliaries[i].name+' is a vertical integral, but it is needed for the calculation of another vertical integral. This does not work.'+chr(13);
  //second consistency check: No auxiliary variable with vertLoc=WAT or SUR and myCalcBeforeZIntegral=1 may depend on an auxiliary variable with vertLoc=SED
  for i:=length(auxiliaries)-1 downto 0 do
    if (auxiliaries[i].vertLoc=0) and (auxiliaries[i].myCalcBeforeZIntegral=1) then
    begin
      for j:=0 to i-1 do
        if (auxiliaries[j].vertLoc=1) or (auxiliaries[j].vertLoc=2) then
        begin
          if  ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[1])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[3])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[4])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[5])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[6])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[7])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[8])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[9])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2])
           or ContainedInString(auxiliaries[j].name,auxiliaries[i].formula) then
           s:=s+'Auxiliary variable '+auxiliaries[i].name+' has vertLoc=WAT and depends on auxiliary variable '+auxiliaries[j].name+'with vertLoc=SED or SUR. This is normally fine because the latter one can be precalculated, but in this case, a precalculation of the former one is also needed. More than one additional k-loop for precalculation is not supported.'+chr(13);
        end;
    end;

  //seek for maximum number of iterations in auxiliaries
  maxIterations :=0;
  for i:=0 to length(auxiliaries)-1 do
    if auxiliaries[i].iterations > maxIterations then
      maxIterations :=auxiliaries[i].iterations;

  //check which processes may take place in the pore water and which shall not:
  //processes take place in the pore water if all of the tracers they consume/produce have isInPorewater=1 or vertLoc=SED
  for i:=0 to length(processes)-1 do
  begin
    processes[i].myIsInPorewater:=1;
    for j:=0 to length(processes[i].input)-1 do
      if processes[i].input[j].myTracerNum>=0 then
        if (tracers[processes[i].input[j].myTracerNum].isInPorewater=0) and (tracers[processes[i].input[j].myTracerNum].isActive=1) and (tracers[processes[i].input[j].myTracerNum].vertLoc<>1) then processes[i].myIsInPorewater:=0;
    for j:=0 to length(processes[i].output)-1 do
      if processes[i].output[j].myTracerNum>=0 then
        if (tracers[processes[i].output[j].myTracerNum].isInPorewater=0) and (tracers[processes[i].output[j].myTracerNum].isActive=1) and (tracers[processes[i].output[j].myTracerNum].vertLoc<>1) then processes[i].myIsInPorewater:=0;
  end;

  //If we have some processes with myIsSelected=true, we need to find out which auxiliaries, constants and limitations are needed to calculate them.
  //These shall then get myIsSelected=true
  //first initialize all with false
  for i:=0 to length(limitations)-1 do
    limitations[i].myIsSelected:=false;
  for i:=0 to length(auxiliaries)-1 do
    auxiliaries[i].myIsSelected:=false;
  for i:=0 to length(constants)-1 do
    constants[i].myIsSelected:=false;
  for i:=0 to length(tracers)-1 do
    tracers[i].myIsSelected:=false;
  ii:=0;
  //now check if we have some processes with myIsSelected=true
  for i:=0 to length(processes)-1 do
    if (processes[i].myIsSelected) and (processes[i].isActive=1) then ii:=ii+1;
  if ii>0 then
  begin
    //first step: All limitations which limit these processes get myIsSelected=true
    for i:=0 to length(processes)-1 do
      if (processes[i].myIsSelected) and (processes[i].isActive=1) then
        for k:=0 to length(processes[i].limitations)-1 do
          if processes[i].limitations[k].limitationNum>=0 then
            limitations[processes[i].limitations[k].limitationNum].myIsSelected:=true;
    //second step: All auxiliaries on which the myIsSelected=true processes depend get myIsSelected=true
    for i:=0 to length(processes)-1 do
      if (processes[i].myIsSelected) and (processes[i].isActive=1) then
        for j:=0 to length(auxiliaries)-1 do
        begin
          if ContainedInString(auxiliaries[j].name,processes[i].turnover) then auxiliaries[j].myIsSelected:=true;
          for k:=0 to length(processes[i].limitations)-1 do
            if ContainedInString(auxiliaries[j].name,limitations[processes[i].limitations[k].limitationNum].value) then auxiliaries[j].myIsSelected:=true;
          for k:=0 to length(processes[i].input)-1 do
            if ContainedInString(auxiliaries[j].name,processes[i].input[k].amount) then auxiliaries[j].myIsSelected:=true;
          for k:=0 to length(processes[i].output)-1 do
            if ContainedInString(auxiliaries[j].name,processes[i].output[k].amount) then auxiliaries[j].myIsSelected:=true;
        end;
    //third step: The same for the constants on which the processes depend
    for i:=0 to length(processes)-1 do
      if (processes[i].myIsSelected) and (processes[i].isActive=1) then
        for j:=0 to length(constants)-1 do
        begin
          if ContainedInString(constants[j].name,processes[i].turnover) then constants[j].myIsSelected:=true;
          for k:=0 to length(processes[i].limitations)-1 do
            if ContainedInString(constants[j].name,limitations[processes[i].limitations[k].limitationNum].value) then constants[j].myIsSelected:=true;
          for k:=0 to length(processes[i].input)-1 do
            if ContainedInString(constants[j].name,processes[i].input[k].amount) then constants[j].myIsSelected:=true;
          for k:=0 to length(processes[i].output)-1 do
            if ContainedInString(constants[j].name,processes[i].output[k].amount) then constants[j].myIsSelected:=true;
        end;
    //fourth step: The same for the tracers on which the processes depend, or which are changed by them
    for i:=0 to length(processes)-1 do
      if (processes[i].myIsSelected) and (processes[i].isActive=1) then
        for j:=0 to length(tracers)-1 do
        begin
          if ContainedInString(tracers[j].name,processes[i].turnover) then tracers[j].myIsSelected:=true;
          for k:=0 to length(processes[i].limitations)-1 do
            if ContainedInString(tracers[j].name,limitations[processes[i].limitations[k].limitationNum].value) then tracers[j].myIsSelected:=true;
          for k:=0 to length(processes[i].input)-1 do
            if processes[i].input[k].myTracerNum=j then tracers[j].myIsSelected:=true;
          for k:=0 to length(processes[i].output)-1 do
            if processes[i].output[k].myTracerNum=j then tracers[j].myIsSelected:=true;
        end;
    //fifth step: Find auxiliaries on which the other auxiliaries depend
    //this step has to be repeated until no new ones are found, this is saved in ii
    ii:=1;
    while ii>0 do
    begin
      ii:=0;
      for i:=0 to length(auxiliaries)-1 do
        if auxiliaries[i].myIsSelected then
          for j:=0 to length(auxiliaries)-1 do
            if auxiliaries[j].myIsSelected=false then
            begin
              if ContainedInString(auxiliaries[j].name,auxiliaries[i].formula) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[1]) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2]) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[3]) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[4]) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[5]) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[6]) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[7]) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[8]) or
                 ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[9]) then //auxiliary i is dependent on auxiliary j
              begin
                ii:=ii+1;
                auxiliaries[j].myIsSelected:=true;
              end;
            end;
    end;
    //sixth step: Find constants which are needed to calculate the auxiliaries
    for i:=0 to length(auxiliaries)-1 do
      if auxiliaries[i].myIsSelected then
        for j:=0 to length(constants)-1 do
        begin
          if ContainedInString(constants[j].name,auxiliaries[i].formula) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[1]) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[2]) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[3]) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[4]) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[5]) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[6]) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[7]) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[8]) or
             ContainedInString(constants[j].name,auxiliaries[i].temp[9]) then //auxiliary i is dependent on constant j
            constants[j].myIsSelected:=true;
        end;
    //final step: Find tracers which are needed to calculate the auxiliaries
    for i:=0 to length(auxiliaries)-1 do
      if auxiliaries[i].myIsSelected then
        for j:=0 to length(tracers)-1 do
        begin
          if ContainedInString(tracers[j].name,auxiliaries[i].formula) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[1]) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[2]) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[3]) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[4]) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[5]) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[6]) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[7]) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[8]) or
             ContainedInString(tracers[j].name,auxiliaries[i].temp[9]) then //auxiliary i is dependent on constant j
            tracers[j].myIsSelected:=true;
        end;
  end;

  result:=s;
end;

function CheckAuxiliaryOrder: String;
var
  i,j: Integer;
  s: String;
begin
  s:='';
  //perform a consistency check: auxiliaries need to have dependencies in the correct order,
  //unless a) the latter has useElsewhere=1 or
  //       b) both have iterations > 0
  for i:=0 to length(auxiliaries)-1 do
    for j:=i+1 to length(auxiliaries)-1 do
    begin
      if ContainedInString(auxiliaries[j].name,auxiliaries[i].formula) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[1]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[3]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[4]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[5]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[6]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[7]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[8]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[9]) then //auxiliary i is dependent on auxiliary j
      begin
        if (auxiliaries[j].iterations=0) then  //the required variable is not iterative (and therefore not pre-calculated in the iterative loop)
          if auxiliaries[j].isUsedElsewhere=0 then // the value does not come from the last time step
            s:=s+'Auxiliary variable '+auxiliaries[i].name+' depends on the value of auxiliary variable '+auxiliaries[j].name+', but they are calculated in the wrong order.'+chr(13);
      end;
    end;

  //perform another consistency check: auxiliaries and processes which have vertLoc=WAT may not depend on those with vertLoc=SED or vertLoc=SUR
  //unless the latter has isZIntegral=1
  for i:=0 to length(auxiliaries)-1 do
    for j:=0 to length(auxiliaries)-1 do
    begin
      if ContainedInString(auxiliaries[j].name,auxiliaries[i].formula) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[1]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[3]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[4]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[5]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[6]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[7]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[8]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[9]) then //auxiliary i is dependent on auxiliary j
      begin
        if (auxiliaries[j].vertLoc=1) or (auxiliaries[j].vertLoc=2) then //the required variable has vertLoc=SED or vertLoc=SUR
          if auxiliaries[i].vertLoc=0 then //the requiring variable has vertLoc=WAT
            if auxiliaries[j].isZIntegral<>1 then
              s:=s+'Auxiliary variable '+auxiliaries[i].name+' with vertLoc=WAT depends on the value of auxiliary variable '+auxiliaries[j].name+', which has vertLoc=SED or vertLoc=SUR.'+chr(13);
      end;
    end;

  //perform another consistency check: auxiliaries and processes which have vertLoc=WAT may not depend on those with vertLoc=SED or vertLoc=SUR
  //unless a) the latter has useElsewhere=1 or
  //       b) both have iterations > 0
  for i:=0 to length(processes)-1 do
    for j:=0 to length(auxiliaries)-1 do
    begin
      if ContainedInString(auxiliaries[j].name,processes[i].turnover) then //auxiliary i is dependent on auxiliary j
      begin
        if (auxiliaries[j].vertLoc=1) or (auxiliaries[j].vertLoc=2) then //the required variable has vertLoc=SED or vertLoc=SUR
          if Processes[i].vertLoc=0 then //the requiring variable has vertLoc=WAT
            s:=s+'Process '+processes[i].name+' with vertLoc=WAT depends on the value of auxiliary variable '+auxiliaries[j].name+', which has vertLoc=SED or vertLoc=SUR.'+chr(13);
      end;
    end;
  result:=s;
end;

function ElementCreatedByProcess_sink(element: Integer; process: Integer; ConsiderVirtualTracers: Boolean=true):Real;
var
  sum: Real;
  i, j, k, l: Integer;
begin
  sum:=0;
  i:=process;
  for j:=0 to length(processes[i].input)-1 do
  begin
    k:=processes[i].input[j].myTracerNum;
    if (tracers[k].isActive=1) or ConsiderVirtualTracers then
      for l:=0 to length(tracers[k].contents)-1 do
        if tracers[k].contents[l].myElementNum = element then
          sum:=sum+tracers[k].contents[l].myAmount*processes[i].input[j].myAmount;
  end;
  result:=sum;
end;

function ElementCreatedByProcess_source(element: Integer; process: Integer; ConsiderVirtualTracers: Boolean=true):Real;
var
  sum: Real;
  i, j, k, l: Integer;
begin
  sum:=0;
  i:=process;
  for j:=0 to length(processes[i].output)-1 do
  begin
    k:=processes[i].output[j].myTracerNum;
    if (tracers[k].isActive=1) or ConsiderVirtualTracers then
      for l:=0 to length(tracers[k].contents)-1 do
        if tracers[k].contents[l].myElementNum = element then
          sum:=sum+tracers[k].contents[l].myAmount*processes[i].output[j].myAmount;
  end;
  result:=sum;
end;

function ElementCreatedByProcess(element: Integer; process: Integer; ConsiderVirtualTracers: Boolean=true):Real;
var
  sum: Real;
begin
  sum := ElementCreatedByProcess_source(element,process,considerVirtualTracers) - ElementCreatedByProcess_sink(element,process,considerVirtualTracers);
  result:=sum;
end;

function ProcessIsConservative(process: Integer; ConsiderVirtualTracers: Boolean=true):Boolean;
const
  epsilon=0.00000000001;
var
  i: Integer;
  sum: Real;
begin
  result:=true;
  for i:=0 to length(elements)-1 do
  begin
    sum:=ElementCreatedByProcess(i,process,ConsiderVirtualTracers);
    if sum > epsilon then
      result:=false
    else if sum < -epsilon then
      result:=false;
  end;
end;

procedure FindSourcesSinks(element: Integer; var sources: TStringList; var sinks: TStringList; ConsiderVirtualTracers: Boolean=true);
const
  epsilon=0.00000000001;
var
  i: Integer;
  sum: real;
begin
  sources.Clear;
  sinks.Clear;
  for i:=0 to length(processes)-1 do
  begin
    sum:=ElementCreatedByProcess(element,i,ConsiderVirtualTracers);
    if sum > epsilon then
      sources.Add(processes[i].name+' ('+processes[i].description+')')
    else if sum < -epsilon then
      sinks.Add(processes[i].name+' ('+processes[i].description+')');
  end;
end;

procedure AutoLimitProcesses;
//Adds the limitation "HARD tracer > 0" to all processes that consume "tracer"
//and have no limitation on "tracer" yet.
//Also, "tracer" has to have set isPositive = 1
//                               isActive   = 1
var
  p, l, i, t: Integer;
  found: Boolean;
  prolim: TErgProcessLimitation;
begin
  for p:=0 to length(processes)-1 do
  begin
    for i:=0 to length(processes[p].input)-1 do
    begin
      t:=processes[p].input[i].myTracerNum;
      if (tracers[t].isPositive=1)
         and (tracers[t].isActive=1)
         and (tracers[t].vertLoc<>3) then
      begin
        //this tracer is a precursor of the process p
        //now, search for existing limitations
        found:=false;
        for l:=0 to length(processes[p].limitations)-1 do
          if lowercase(limitations[processes[p].limitations[l].limitationNum].tracer) = lowercase(tracers[t].name) then
            if processes[p].limitations[l].tracerIsSmall=0 then //process is switched of at small tracer concentrations, not at large ones
              found:=true;
        if not found then //no limitation existed yet, so create the hard one
        begin
          if GetProcessLimitation('HARD '+tracers[t].name+' > 0.0',prolim) then
          begin
            setLength(Processes[p].limitations,length(Processes[p].limitations)+1);
            Processes[p].limitations[length(Processes[p].limitations)-1]:=prolim;
          end;
        end;
      end;
    end;
  end;
  GenerateIndexes;
end;

procedure ApplyLimitations;
//Modifies the process rates by multiplying them with the values of the limitation functions,
// rate := lim_tracer_XXX * old_rate
//Replacement processes get the remaining rate, (1-lim_tracer_XXX) * old_rate.
var
  i,j: Integer;
begin
  for i:=0 to length(processes)-1 do
  begin
    if length(processes[i].limitations)>0 then
      processes[i].turnover:='('+processes[i].turnover+')';
    for j:=0 to length(processes[i].limitations)-1 do
    begin
      if processes[i].limitations[j].tracerIsSmall=0 then
      begin
        if processes[i].limitations[j].myElseProcessNum <> -1 then
          processes[processes[i].limitations[j].myElseProcessNum].turnover :=
            '('+processes[processes[i].limitations[j].myElseProcessNum].turnover+')+'+
            processes[i].limitations[j].myElseProcessRatio+'*(1.0-lim_'+limitations[processes[i].limitations[j].limitationNum].tracer+'_'+IntToStr(processes[i].limitations[j].limitationNum)+')*'+processes[i].turnover;
        processes[i].turnover := processes[i].turnover + '*lim_'+limitations[processes[i].limitations[j].limitationNum].tracer+'_'+IntToStr(processes[i].limitations[j].limitationNum);
      end
      else
      begin
        if processes[i].limitations[j].myElseProcessNum <> -1 then
          processes[processes[i].limitations[j].myElseProcessNum].turnover :=
            '('+processes[processes[i].limitations[j].myElseProcessNum].turnover+')+'+
            processes[i].limitations[j].myElseProcessRatio+'*(lim_'+limitations[processes[i].limitations[j].limitationNum].tracer+'_'+IntToStr(processes[i].limitations[j].limitationNum)+')*'+processes[i].turnover;
        processes[i].turnover := processes[i].turnover + '*(1.0-lim_'+limitations[processes[i].limitations[j].limitationNum].tracer+'_'+IntToStr(processes[i].limitations[j].limitationNum)+')';
      end;
    end;
  end;
end;

procedure GetProcessTypes(var sl: TStringList);
//Fills the string list sl as follows:
//For all occuring process types (processes[i].processType), it contains an entry of this name.
//If it is disabled in modelinfos.inactiveProcessTypes, it is preceded by "0", else by "1".
//e.g. 0gas_exchange
//     1standard
var
  i, p: Integer;
  found: Boolean;
  s, s1: String;
begin
  //first, store all process types in sl, preceded by "1".
  sl.Clear;
  for p:=0 to length(processes)-1 do
  begin
    found:=false;    //make sure each one is listed only once
    for i:=0 to sl.Count-1 do
      if trim(lowercase(processes[p].processType)) = lowercase(copy(sl[i],2,length(sl[i]))) then found:=true;
    if not found then
      sl.Add('1'+trim(processes[p].processType));
  end;
  //now, replace 1 by 0 if processType occurs in modelinfos.inactiveProcessTypes
  s:=modelinfos.inactiveProcessTypes;
  while not (s='') do
  begin
    s1:=trim(semiItem(s));
    for i:=0 to sl.Count-1 do
      if s1 = lowercase(copy(sl[i],2,length(sl[i]))) then
        sl[i]:='0'+copy(sl[i],2,length(sl[i]));
  end;
end;

procedure SetDisabledProcessesInactive(sl: TStringList);
var
  i, p: Integer;
begin
  for p:=0 to length(processes)-1 do
  begin
    for i:=0 to sl.Count-1 do
      if trim(lowercase(processes[p].processType)) = lowercase(copy(sl[i],2,length(sl[i]))) then
      begin
        if copy(sl[i],1,1)='0' then
          processes[p].isActive:=0;
      end;
  end;
end;

procedure SortAuxiliaries;
var
  i,j: Integer;
  L,S: TStringList;
  dependent: Array of Array of Boolean;
  isReady: Boolean;
  sorted_auxiliaries: Array of TErgAuxiliary;
begin
  //Reorder auxiliary variables such that they are calculated in the correct order (their dependencies are resolved).
  //auxiliaries need to have dependencies in the correct order,
  //unless a) the latter has useElsewhere=1 or
  //       b) both have iterations > 0

  //first, create dependency matrix
  setLength(dependent,length(auxiliaries),length(auxiliaries));
  for i:=0 to length(auxiliaries)-1 do
    for j:=0 to length(auxiliaries)-1 do
    begin
      dependent[i,j]:=false;
      if ContainedInString(auxiliaries[j].name,auxiliaries[i].formula) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[1]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[2]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[3]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[4]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[5]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[6]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[7]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[8]) or
         ContainedInString(auxiliaries[j].name,auxiliaries[i].temp[9]) then //auxiliary i is dependent on auxiliary j
        if (auxiliaries[j].iterations=0) then  //the required variable is not iterative (and therefore not pre-calculated in the iterative loop)
          if auxiliaries[j].isUsedElsewhere=0 then // the value does not come from the last time step
            dependent[i,j]:=true;
    end;

  //second, check if a re-sorting is required. This is the case if a former auxiliary is dependent on a latter one.
  isReady:=true;
  for i:=0 to length(auxiliaries)-2 do
    for j:=i+1 to length(auxiliaries)-1 do
      if dependent[i,j] then isReady:=false;

  //third, if necessary, use the Kahn 1962 algorithm to get them sorted
  //Kahn, Arthur B. (1962), "Topological sorting of large networks", Communications of the ACM 5 (11): 558–562, doi:10.1145/368996.369025
  if isReady=false then
  begin
    L:=TStringList.Create;
    S:=TStringList.Create;

    for i:=0 to length(auxiliaries)-1 do
    begin
      isReady:=true;
      for j:=0 to length(auxiliaries)-1 do
        if dependent[i,j] then isReady:=false;
      if isReady then S.Add(IntToStr(i));
    end;

    while S.Count>0 do
    begin
      L.Add(s[0]);
      for i:=0 to length(auxiliaries)-1 do
        if dependent[i,StrToInt(S[0])] then
        begin
          dependent[i,StrToInt(S[0])]:=false;
          isReady:=true;
          for j:=0 to length(auxiliaries)-1 do
            if dependent[i,j] then isReady:=false;
          if isReady then S.Add(IntToStr(i));
        end;
      S.Delete(0);
    end;

    if L.Count < length(auxiliaries) then
      ShowMessage('An error occurred: The auxiliary variables have cyclic dependencies. They were therefore not sorted in the correct order as planned.')
    else
    begin
      setLength(sorted_auxiliaries,length(auxiliaries));        //temporary list for re-sorting
      for i:=0 to length(auxiliaries)-1 do
        copyErgAuxiliary(auxiliaries[StrToInt(L[i])],sorted_auxiliaries[i]);
      for i:=0 to length(auxiliaries)-1 do
        copyErgAuxiliary(sorted_auxiliaries[i],auxiliaries[i]);
      ShowMessage('Auxiliary variables were succesfully sorted such that all of them are calculated in the correct order.')
    end;

    L.Free; S.Free;
  end;
end;

function FinishLoadingAddOn: String;
var
  i,j: Integer;
  found: Boolean;
  delete_name: String;
  sorted_auxiliaries: Array of TErgAuxiliary;
  dependent: Array of Array of Boolean;
  L,S: TStringList;
  isReady: Boolean;
begin
  result:='';
  //first, delete all constants, processes, ... which have to be deleted for the add-on

  found:=true;
  while found do //repeat until nothing more to delete
  begin
    found:=false;
    for i:=0 to length(constants)-1 do                          // seek for a command to delete a constant
      if copy(constants[i].name,1,17)='delete_in_add_on_' then  // found one
      begin
        found:=true;
        delete_name:=trim(lowercase(copy(constants[i].name,18,length(constants[i].name))));
        for j:=i to length(constants)-2 do                      // delete this command from the list of constants
          copyErgConstant(constants[j+1],constants[j]);
        setLength(constants,length(constants)-1);
        break;
      end;
    if found then                                               // now seek the constant
    begin
      i:=0;
      while (trim(lowercase(constants[i].name)) <> delete_name) and (i<length(constants)) do
        i:=i+1;
      if i<length(constants) then                               // constant has been found
      begin
        for j:=i to length(constants)-2 do                      // delete the constant from the list
          copyErgConstant(constants[j+1],constants[j]);
        setLength(constants,length(constants)-1);
      end
      else
        result:=result+'Constant '+delete_name+' should be deleted by the add-on, but was not found. Ignored.'+chr(13);
    end;
  end;

  found:=true;
  while found do //repeat until nothing more to delete
  begin
    found:=false;
    for i:=0 to length(Tracers)-1 do                          // seek for a command to delete a Tracer
      if copy(Tracers[i].name,1,17)='delete_in_add_on_' then  // found one
      begin
        found:=true;
        delete_name:=trim(lowercase(copy(Tracers[i].name,18,length(Tracers[i].name))));
        for j:=i to length(Tracers)-2 do                      // delete this command from the list of Tracers
          copyErgTracer(Tracers[j+1],Tracers[j]);
        setLength(Tracers,length(Tracers)-1);
        break;
      end;
    if found then                                               // now seek the Tracer
    begin
      i:=0;
      while (trim(lowercase(Tracers[i].name)) <> delete_name) and (i<length(Tracers)) do
        i:=i+1;
      if i<length(Tracers) then                               // Tracer has been found
      begin
        for j:=i to length(Tracers)-2 do                      // delete the Tracer from the list
          copyErgTracer(Tracers[j+1],Tracers[j]);
        setLength(Tracers,length(Tracers)-1);
      end
      else
        result:=result+'Tracer '+delete_name+' should be deleted by the add-on, but was not found. Ignored.'+chr(13);
    end;
  end;

  found:=true;
  while found do //repeat until nothing more to delete
  begin
    found:=false;
    for i:=0 to length(Auxiliaries)-1 do                          // seek for a command to delete a Auxiliary
      if copy(Auxiliaries[i].name,1,17)='delete_in_add_on_' then  // found one
      begin
        found:=true;
        delete_name:=trim(lowercase(copy(Auxiliaries[i].name,18,length(Auxiliaries[i].name))));
        for j:=i to length(Auxiliaries)-2 do                      // delete this command from the list of Auxiliaries
          copyErgAuxiliary(Auxiliaries[j+1],Auxiliaries[j]);
        setLength(Auxiliaries,length(Auxiliaries)-1);
        break;
      end;
    if found then                                               // now seek the Auxiliary
    begin
      i:=0;
      while (trim(lowercase(Auxiliaries[i].name)) <> delete_name) and (i<length(Auxiliaries)) do
        i:=i+1;
      if i<length(Auxiliaries) then                               // Auxiliary has been found
      begin
        for j:=i to length(Auxiliaries)-2 do                      // delete the Auxiliary from the list
          copyErgAuxiliary(Auxiliaries[j+1],Auxiliaries[j]);
        setLength(Auxiliaries,length(Auxiliaries)-1);
      end
      else
        result:=result+'Auxiliary '+delete_name+' should be deleted by the add-on, but was not found. Ignored.'+chr(13);
    end;
  end;

  found:=true;
  while found do //repeat until nothing more to delete
  begin
    found:=false;
    for i:=0 to length(Processes)-1 do                          // seek for a command to delete a Process
      if copy(Processes[i].name,1,17)='delete_in_add_on_' then  // found one
      begin
        found:=true;
        delete_name:=trim(lowercase(copy(Processes[i].name,18,length(Processes[i].name))));
        for j:=i to length(Processes)-2 do                      // delete this command from the list of Processes
          copyErgProcess(Processes[j+1],Processes[j]);
        setLength(Processes,length(Processes)-1);
        break;
      end;
    if found then                                               // now seek the Process
    begin
      i:=0;
      while (trim(lowercase(Processes[i].name)) <> delete_name) and (i<length(Processes)) do
        i:=i+1;
      if i<length(Processes) then                               // Process has been found
      begin
        for j:=i to length(Processes)-2 do                      // delete the Process from the list
          copyErgProcess(Processes[j+1],Processes[j]);
        setLength(Processes,length(Processes)-1);
      end
      else
        result:=result+'Process '+delete_name+' should be deleted by the add-on, but was not found. Ignored.'+chr(13);
    end;
  end;

  found:=true;
  while found do //repeat until nothing more to delete
  begin
    found:=false;
    for i:=0 to length(Elements)-1 do                          // seek for a command to delete a Element
      if copy(Elements[i].name,1,17)='delete_in_add_on_' then  // found one
      begin
        found:=true;
        delete_name:=trim(lowercase(copy(Elements[i].name,18,length(Elements[i].name))));
        for j:=i to length(Elements)-2 do                      // delete this command from the list of Elements
          copyErgElement(Elements[j+1],Elements[j]);
        setLength(Elements,length(Elements)-1);
        break;
      end;
    if found then                                               // now seek the Element
    begin
      i:=0;
      while (trim(lowercase(Elements[i].name)) <> delete_name) and (i<length(Elements)) do
        i:=i+1;
      if i<length(Elements) then                               // Element has been found
      begin
        for j:=i to length(Elements)-2 do                      // delete the Element from the list
          copyErgElement(Elements[j+1],Elements[j]);
        setLength(Elements,length(Elements)-1);
      end
      else
        result:=result+'Element '+delete_name+' should be deleted by the add-on, but was not found. Ignored.'+chr(13);
    end;
  end;

  found:=true;
  while found do //repeat until nothing more to delete
  begin
    found:=false;
    for i:=0 to length(CElements)-1 do                          // seek for a command to delete a CElement
      if copy(CElements[i].color,1,17)='delete_in_add_on_' then  // found one
      begin
        found:=true;
        delete_name:=trim(lowercase(copy(CElements[i].color,18,length(CElements[i].color))))+'_'+trim(lowercase(celements[i].element));
        for j:=i to length(CElements)-2 do                      // delete this command from the list of CElements
          copyErgCElement(CElements[j+1],CElements[j]);
        setLength(CElements,length(CElements)-1);
        break;
      end;
    if found then                                               // now seek the CElement
    begin
      i:=0;
      while (trim(lowercase(CElements[i].color))+'_'+trim(lowercase(CElements[i].element)) <> delete_name) and (i<length(CElements)) do
        i:=i+1;
      if i<length(CElements) then                               // CElement has been found
      begin
        for j:=i to length(CElements)-2 do                      // delete the CElement from the list
          copyErgCElement(CElements[j+1],CElements[j]);
        setLength(CElements,length(CElements)-1);
      end
      else
        result:=result+'CElement '+delete_name+' should be deleted by the add-on, but was not found. Ignored.'+chr(13);
    end;
  end;

  sortAuxiliaries;
end;

end.
