unit erg_base;

interface

uses erg_types, classes, Dialogs;

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;
  cellheightTimesDensity: String;
  maxIterations: Integer;

function LoadConstants(filename: String):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);

function LoadElements(filename: String):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);

function LoadTracers(filename: String):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);

function LoadAuxiliaries(filename: String):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);

function LoadProcesses(filename: String):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);

function LoadCElements(filename: String):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);

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: String;   //fills all the "my..." values in the records
                                    //returns error string or ''

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(i: Integer):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 initKeywords;

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,25); setLength(keywordDescriptions,25);
  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]';
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;

//************************** 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);
  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
          constants[length(constants)-1].name:=s2
        else if lowercase(s1)='value' then
        begin
          constants[length(constants)-1].valueString:=s2;
          //if pos(';',s2)<=0 then constants[length(constants)-1].value:=StrToFloat(s2);
        end
        else if lowercase(s1)='description' then
          constants[length(constants)-1].description:=s2
        else if lowercase(s1)='comment' then
          constants[length(constants)-1].comment:=s2
        else result:=2; //undefined variables were set
      end;
  end;    //while ended, now check if a constant with this name exists
  if constants[length(constants)-1].name = '' then
    setLength(constants,length(constants)-1)
  else
  begin
    p:=-1;
    //seek the constant with the same name
    for i:=0 to length(constants)-2 do
      if lowercase(constants[i].name) = lowercase(constants[length(constants)-1].name) then
        p:=i;
    if p>=0 then  //process found
    begin
      //apply the settings of the loaded constant to the existing one
      copyErgConstant(constants[length(constants)-1],constants[p]);
      //remove the loaded process
      setLength(constants,length(constants)-1);
    end;
  end;
end;

function LoadConstants(filename: String):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DecimalSeparator:='.';
  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,'!   comment=        comment, e.g. how certain this value is, literature,..., default=""');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleConstant(var F: TextFile; i: Integer);
var c: TErgConstant;
begin
  InitErgConstant(c);
      writeln(F,'name        = '+constants[i].name);
      writeln(F,'value       = '+constants[i].valueString);
      if constants[i].description <> c.description then
        writeln(F,'description = '+constants[i].description);
      if constants[i].comment <> c.comment then
        writeln(F,'comment     = '+constants[i].comment);
end;

procedure SaveConstants(filename: String);
var F: TextFile;
    i: Integer;
begin
  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);
  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
          Elements[length(Elements)-1].name:=s2
        else if lowercase(s1)='description' then
          Elements[length(Elements)-1].description:=s2
        else if lowercase(s1)='comment' then
          Elements[length(Elements)-1].comment:=s2
        else result:=2; //undefined variables were set
      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)
  else
  begin
    p:=-1;
    //seek the Element with the same name
    for i:=0 to length(Elements)-2 do
      if lowercase(Elements[i].name) = lowercase(Elements[length(Elements)-1].name) then
        p:=i;
    if p>=0 then  //process found
    begin
      //apply the settings of the loaded Element to the existing one
      copyErgElement(Elements[length(Elements)-1],Elements[p]);
      //remove the loaded process
      setLength(Elements,length(Elements)-1);
    end;
  end;
end;

function LoadElements(filename: String):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DecimalSeparator:='.';
  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 SaveSingleElement(var F: TextFile; i: Integer);
var c: TErgElement;
begin
  InitErgElement(c);
      writeln(F,'name        = '+elements[i].name);
      writeln(F,'description = '+elements[i].description);
      if elements[i].comment<>c.comment then
        writeln(F,'comment = '+elements[i].comment);
end;

procedure SaveElements(filename: String);
var F: TextFile;
    i: Integer;
begin
  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, p: Integer;
begin
  result:=0;
  setLength(Tracers,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
          Tracers[length(Tracers)-1].name:=s2
        else if lowercase(s1)='description' then
          Tracers[length(Tracers)-1].description:=s2
        else if lowercase(s1)='comment' then
          Tracers[length(Tracers)-1].comment:=s2
        else if lowercase(s1)='solubility' then
          Tracers[length(Tracers)-1].solubility:=s2
        else if lowercase(s1)='schmidtnumber' then
          Tracers[length(Tracers)-1].schmidtNumber:=s2
        else if lowercase(s1)='gasname' then
          Tracers[length(Tracers)-1].gasName:=s2
        else if lowercase(s1)='molarmass' then
          Tracers[length(Tracers)-1].molarMass:=s2
        else if lowercase(s1)='contents' then
        begin
          SetLength(Tracers[length(Tracers)-1].contents,StrToInt(s2));
          for i:=0 to length(Tracers[length(Tracers)-1].contents)-1 do
          begin
            s1:='';
            while s1='' do
            begin
              readln(F,s);
              splitEquation(s,s1,s2);
            end;
            Tracers[length(Tracers)-1].contents[i].element:=s1;
            Tracers[length(Tracers)-1].contents[i].amount:=s2;
          end;
        end
        else if lowercase(s1)='verticaldistribution' then
          Tracers[length(Tracers)-1].verticalDistribution:=s2
        else if lowercase(s1)='vertspeed' then
          Tracers[length(Tracers)-1].vertSpeed:=s2
        else if lowercase(s1)='vertdiff' then
          Tracers[length(Tracers)-1].vertDiff:=s2
        else if lowercase(s1)='opacity' then
          Tracers[length(Tracers)-1].opacity:=s2
        else if (lowercase(s1)='isflat') or (lowercase(s1)='vertloc') then
          Tracers[length(Tracers)-1].vertLoc:=StrToIntVertLoc(s2)
        else if lowercase(s1)='ispositive' then
          Tracers[length(Tracers)-1].isPositive:=StrToInt(s2)
        else if lowercase(s1)='ismixed' then
          Tracers[length(Tracers)-1].isMixed:=StrToInt(s2)
        else if lowercase(s1)='iscombined' then
          Tracers[length(Tracers)-1].isCombined:=StrToInt(s2)
        else if lowercase(s1)='isstiff' then
          Tracers[length(Tracers)-1].isStiff:=StrToInt(s2)
        else if lowercase(s1)='initvalue' then
          Tracers[length(Tracers)-1].initValue:=s2
        else if lowercase(s1)='useinitvalue' then
          Tracers[length(Tracers)-1].useInitValue:=StrToInt(s2)
        else if lowercase(s1)='childof' then
          Tracers[length(Tracers)-1].childOf:=s2
        else if lowercase(s1)='dimension' then
          Tracers[length(Tracers)-1].dimension:=StrToInt(s2)
        else if lowercase(s1)='masslimits' then
          Tracers[length(Tracers)-1].massLimits:=s2
        else if lowercase(s1)='atmosdep' then
          Tracers[length(Tracers)-1].atmosDep:=StrToInt(s2)
        else if lowercase(s1)='riverdep' then
          Tracers[length(Tracers)-1].riverDep:=StrToInt(s2)
        else if lowercase(s1)='isoutput' then
          Tracers[length(Tracers)-1].isOutput:=StrToInt(s2)
        else if lowercase(s1)='isactive' then
          Tracers[length(Tracers)-1].isActive:=StrToInt(s2)
        else
          result:=2; //undefined variables were set
      end;
  end;    //while ended, now check if a Tracer with this name exists
  if Tracers[length(Tracers)-1].name = '' then
    setLength(Tracers,length(Tracers)-1)
  else
  begin
    p:=-1;
    //seek the Tracer with the same name
    for i:=0 to length(Tracers)-2 do
      if lowercase(Tracers[i].name) = lowercase(Tracers[length(Tracers)-1].name) then
        p:=i;
    if p>=0 then  //process found
    begin
      //apply the settings of the loaded Tracer to the existing one
      copyErgTracer(Tracers[length(Tracers)-1],Tracers[p]);
      //remove the loaded process
      setLength(Tracers,length(Tracers)-1);
    end;
  end;
end;

function LoadTracers(filename: String):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DecimalSeparator:='.';
  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,'! 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,'! isStiff=        0=not stiff (default); 1=stiff, use Patankar method if concentration declines; 2=stiff, always use modified Patankar method');
  writeln(F,'! isMixed=        1=mix with neighbour cells if negative, 0=do not, default=1, only applicable to tracers with vertLoc=WAT');
  writeln(F,'! isOutput=       1=occurs as output in model results (default); 0=only internal use');
  writeln(F,'! isPositive=     0=may be negative, 1=always positive (default)');
  writeln(F,'! massLimits=     semicolon-seperated string of (dimension-1) FORTRAN expressions for mass limits for stage-resolving models [mmol], default=""');
  writeln(F,'! opacity=        fortran formula for opacity (for light limitation), 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,'! schmidtNumber=  name of an auxiliary variable describing the Schmidt number [1] for gasses which flow through the surface, 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,'! riverDep=       0=no river depositon of this tracer, 1=river deposition of this tracer (default)');
  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, 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,'! comment=        e.g. "represents a certain kind of phytoplankton", default=""');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleTracer(var F: TextFile; i: Integer);
var
  c: TErgTracer;
  j: Integer;
begin
  InitErgTracer(c);
      writeln(F,'name              = '+Tracers[i].name);
      writeln(F,'description       = '+Tracers[i].description);
      if Tracers[i].comment <> c.comment then
        writeln(F,'comment           = '+Tracers[i].comment);
      if length(tracers[i].contents)>0 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].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));
end;

procedure SaveTracers(filename: String);
var F: TextFile;
    i: Integer;
begin
  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);
  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
          auxiliaries[length(auxiliaries)-1].name:=s2
        else if lowercase(s1)='temp1' then
          auxiliaries[length(auxiliaries)-1].temp[1]:=s2
        else if lowercase(s1)='temp2' then
          auxiliaries[length(auxiliaries)-1].temp[2]:=s2
        else if lowercase(s1)='temp3' then
          auxiliaries[length(auxiliaries)-1].temp[3]:=s2
        else if lowercase(s1)='temp4' then
          auxiliaries[length(auxiliaries)-1].temp[4]:=s2
        else if lowercase(s1)='temp5' then
          auxiliaries[length(auxiliaries)-1].temp[5]:=s2
        else if lowercase(s1)='temp6' then
          auxiliaries[length(auxiliaries)-1].temp[6]:=s2
        else if lowercase(s1)='temp7' then
          auxiliaries[length(auxiliaries)-1].temp[7]:=s2
        else if lowercase(s1)='temp8' then
          auxiliaries[length(auxiliaries)-1].temp[8]:=s2
        else if lowercase(s1)='temp9' then
          auxiliaries[length(auxiliaries)-1].temp[9]:=s2
        else if lowercase(s1)='formula' then
          auxiliaries[length(auxiliaries)-1].formula:=s2
        else if lowercase(s1)='calcafterprocesses' then
          auxiliaries[length(auxiliaries)-1].calcAfterProcesses:=strToInt(s2)
        else if lowercase(s1)='iterations' then
          auxiliaries[length(auxiliaries)-1].iterations:=StrToInt(s2)
        else if lowercase(s1)='iterinit' then
          auxiliaries[length(auxiliaries)-1].iterInit:=s2
        else if lowercase(s1)='isusedelsewhere' then
          auxiliaries[length(auxiliaries)-1].isUsedElsewhere:=strToInt(s2)
        else if lowercase(s1)='iszgradient' then
          auxiliaries[length(auxiliaries)-1].isZGradient:=StrToInt(s2)
        else if lowercase(s1)='iszintegral' then
          auxiliaries[length(auxiliaries)-1].isZIntegral:=StrToInt(s2)
        else if lowercase(s1)='description' then
          auxiliaries[length(auxiliaries)-1].description:=s2
        else if lowercase(s1)='comment' then
          auxiliaries[length(auxiliaries)-1].comment:=s2
        else if (lowercase(s1)='isflat') or (lowercase(s1)='vertloc') then
          auxiliaries[length(auxiliaries)-1].vertLoc:=strToIntVertLoc(s2)
        else if lowercase(s1)='isoutput' then
          auxiliaries[length(auxiliaries)-1].isOutput:=strToInt(s2)
        else result:=2; //undefined variables were set
      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)
  else
  begin
    p:=-1;
    //seek the Auxiliary with the same name
    for i:=0 to length(Auxiliaries)-2 do
      if lowercase(Auxiliaries[i].name) = lowercase(Auxiliaries[length(Auxiliaries)-1].name) then
        p:=i;
    if p>=0 then  //process found
    begin
      //apply the settings of the loaded Auxiliary to the existing one
      copyErgAuxiliary(Auxiliaries[length(Auxiliaries)-1],Auxiliaries[p]);
      //remove the loaded process
      setLength(Auxiliaries,length(Auxiliaries)-1);
    end;
  end;
end;

function LoadAuxiliaries(filename: String):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DecimalSeparator:='.';
  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,'!');
  writeln(F,'! All entries with the same value of calcAfterProcesses are calculated in given order.');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleAuxiliary(var F: TextFile; i: Integer);
var c: TErgAuxiliary;
    j: Integer;
begin
  InitErgAuxiliary(c);
      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]);
      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 then
        writeln(F,'description        = '+auxiliaries[i].description);
      if auxiliaries[i].comment <> c.comment then
        writeln(F,'comment            = '+auxiliaries[i].comment);
      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));
end;

procedure SaveAuxiliaries(filename: String);
var F: TextFile;
    i: Integer;
begin
  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(i: Integer):String;
var
  s: String;
  j: Integer;
begin
      s:='';
      for j:=0 to length(Processes[i].input)-1 do
      begin
        if j>0 then s:=s + ' + ';
        if (Processes[i].input[j].amount='1') or (Processes[i].input[j].amount='1.0') then
          s:=s+Processes[i].input[j].tracer
        else
          s:=s+Processes[i].input[j].amount+'*'+Processes[i].input[j].tracer;
      end;
      s:=s+' -> ';
      for j:=0 to length(Processes[i].output)-1 do
      begin
        if j>0 then s:=s + ' + ';
        if (Processes[i].output[j].amount='1') or (Processes[i].output[j].amount='1.0') then
          s:=s+Processes[i].output[j].tracer
        else
          s:=s+Processes[i].output[j].amount+'*'+Processes[i].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,p: Integer;
  prolim:TErgProcessLimitation;
begin
  result:=0;
  setLength(Processes,length(processes)+1);
  initErgProcess(Processes[length(Processes)-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       //Process ended
        break
      else if lowercase(s1)='name' then
        Processes[length(Processes)-1].name:=s2
      else if lowercase(s1)='description' then
        Processes[length(Processes)-1].description:=s2
      else if lowercase(s1)='comment' then
        Processes[length(Processes)-1].comment:=s2
      else if (lowercase(s1)='rate') or (lowercase(s1)='turnover') then
        Processes[length(Processes)-1].turnover:=s2
      else if lowercase(s1)='comment' then
        Processes[length(Processes)-1].comment:=s2
      else if lowercase(s1)='input' then
      begin
        SetLength(Processes[length(Processes)-1].input,StrToInt(s2));
        for i:=0 to length(Processes[length(Processes)-1].input)-1 do
        begin
          s1:='';
          while s1='' do
          begin
            readln(F,s);
            splitEquation(s,s1,s2);
          end;
          Processes[length(Processes)-1].input[i].tracer:=s1;
          Processes[length(Processes)-1].input[i].amount:=s2;
        end;
      end
      else if lowercase(s1)='output' then
      begin
        SetLength(Processes[length(Processes)-1].output,StrToInt(s2));
        for i:=0 to length(Processes[length(Processes)-1].output)-1 do
        begin
          s1:='';
          while s1='' do
          begin
            readln(F,s);
            splitEquation(s,s1,s2);
          end;
          Processes[length(Processes)-1].output[i].tracer:=s1;
          Processes[length(Processes)-1].output[i].amount:=s2;
        end;
      end
      else if lowercase(s1)='equation' then
      begin
        if SetInputOutputEquation(s2,Processes[length(Processes)-1])=false then
          result:=3;
      end
      else if lowercase(s1)='repaint' then
      begin
        SetLength(Processes[length(Processes)-1].repaint,StrToInt(s2));
        for i:=0 to length(Processes[length(Processes)-1].repaint)-1 do
        begin
          s1:='';
          while s1='' do
          begin
            readln(F,s);
            splitEquation(s,s1,s2);
          end;
          splitSpace(s1,s,s1);
          Processes[length(Processes)-1].repaint[i].oldColor:=s;
          Processes[length(Processes)-1].repaint[i].element:=s1;
          Processes[length(Processes)-1].repaint[i].newColor:=s2;
        end;
      end
      else if lowercase(s1)='limitation' then
      begin
        if GetProcessLimitation(s2,prolim) then
        begin
          setLength(Processes[length(Processes)-1].limitations,length(Processes[length(Processes)-1].limitations)+1);
          Processes[length(Processes)-1].limitations[length(Processes[length(Processes)-1].limitations)-1]:=prolim;
        end
        else
          result:=4;
      end
      else if lowercase(s1)='feedingefficiency' then
        Processes[length(Processes)-1].feedingEfficiency:=s2
      else if (lowercase(s1)='isflat') or (lowercase(s1)='vertloc') then
        Processes[length(Processes)-1].vertLoc:=StrToIntVertLoc(s2)
      else if lowercase(s1)='processtype' then
        Processes[length(Processes)-1].processType:=s2
      else if lowercase(s1)='isoutput' then
        Processes[length(Processes)-1].isOutput:=StrToInt(s2)
      else if lowercase(s1)='isactive' then
        Processes[length(Processes)-1].isActive:=StrToInt(s2)
      else result:=2; //undefined variables were set
    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)
  else
  begin
    p:=-1;
    //seek the process with the same name
    for i:=0 to length(processes)-2 do
      if lowercase(processes[i].name) = lowercase(processes[length(processes)-1].name) then
        p:=i;
    if p>=0 then  //process found
    begin
      //apply the settings of the loaded process to the existing one
      copyErgProcess(processes[length(processes)-1],processes[p]);
      //remove the loaded process
      setLength(processes,length(processes)-1);
    end;
  end;
end;

function LoadProcesses(filename: String):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
  DecimalSeparator:='.';
  setLength(limitations,0);
  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,'! 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,'! 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,'! comment=        comment, default=""');
  writeln(F,'!');
  writeln(F,'! Process rates are calculated in the given order.');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleProcess(var F: TextFile; i: Integer);
var
  c: TErgProcess;
  j: Integer;
  s: String;
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(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;
end;

procedure SaveProcesses(filename: String);
var F: TextFile;
    i: Integer;
begin
  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);
  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
          CElements[length(CElements)-1].element:=s2
        else if lowercase(s1)='color' then
          CElements[length(CElements)-1].color:=s2
        else if lowercase(s1)='description' then
          CElements[length(CElements)-1].description:=s2
        else if lowercase(s1)='isaging' then
          CElements[length(CElements)-1].isAging:=s2
        else if lowercase(s1)='atmosdep' then
          CElements[length(CElements)-1].atmosdep:=StrToInt(s2)
        else if lowercase(s1)='riverdep' then
          CElements[length(CElements)-1].riverdep:=StrToInt(s2)
        else if lowercase(s1)='istracer' then
          CElements[length(CElements)-1].isTracer:=s2
        else if lowercase(s1)='comment' then
          CElements[length(CElements)-1].comment:=s2
        else result:=2; //undefined variables were set
      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)
  else
  begin
    p:=-1;
    //seek the CElement with the same element and color
    for i:=0 to length(CElements)-2 do
      if (lowercase(CElements[i].element) = lowercase(CElements[length(CElements)-1].element))
         and (lowercase(CElements[i].color) = lowercase(CElements[length(CElements)-1].color)) then
        p:=i;
    if p>=0 then  //process found
    begin
      //apply the settings of the loaded CElement to the existing one
      copyErgCElement(CElements[length(CElements)-1],CElements[p]);
      //remove the loaded process
      setLength(CElements,length(CElements)-1);
    end;
  end;
end;

function LoadCElements(filename: String):Integer; //returns 0 if succeeded
var F: TextFile;
begin
  DecimalSeparator:='.';
  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,'! 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,'! atmosDep=       1=atmospheric deposition of marked tracers may occur, 0=not (default)');
  writeln(F,'! riverDep=       1=river deposition of marked tracers may occur, 0=not (default)');
  writeln(F,'! comment=        any comments');
  writeln(F,'! *************************************************************************************');
end;

procedure SaveSingleCElement(var F: TextFile; i: Integer);
var c: TErgCElement;
begin
  InitErgCElement(c);
      writeln(F,'element     = '+CElements[i].element);
      writeln(F,'color       = '+CElements[i].color);
      if CElements[i].description <> c.description 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].comment <> c.comment then
        writeln(F,'comment     = '+CElements[i].comment);
end;

procedure SaveCElements(filename: String);
var F: TextFile;
    i: Integer;
begin
  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
  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)='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)='autounixoutput' then
          ModelInfos.autoUnixOutput:=StrToInt(s2)
        else if lowercase(s1)='inactiveprocesstypes' then
          ModelInfos.inactiveProcessTypes:=s2
        else if lowercase(s1)='comment' then
          ModelInfos.comment:=s2
        else result:=2; //undefined variables were set
      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');
  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,'! 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,'! 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
  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.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.autoUnixOutput <> c.autoUnixOutput then
      writeln(F,'autoUnixOutput    = '+IntToStr(ModelInfos.autoUnixOutput));
    if ModelInfos.inactiveProcessTypes <> c.inactiveProcessTypes then
      writeln(F,'inactiveProcessTypes = '+ModelInfos.InactiveProcessTypes);
    if ModelInfos.comment <> c.comment then
      writeln(F,'comment              = '+ModelInfos.comment);
  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
  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: String;
var i,j,k, ii: Integer;
    s: String;
    e: TParseExpression;
    precedingStuff: String;
    numerator, denominator: String;
    ratestring: String;
begin
  s:='';
  DecimalSeparator:='.';

  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].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;

  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;
    for j:=0 to length(tracers)-1 do
      if lowercase(tracers[j].name) = lowercase(tracers[i].childOf) then
        tracers[i].myChildOf:=j;
    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);
    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;
    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
      processes[i].input[j].myTracerNum:=-1;
      for k:=0 to length(tracers)-1 do
        if lowercase(processes[i].input[j].tracer) = lowercase(tracers[k].name) then
          processes[i].input[j].myTracerNum:=k;
      if processes[i].input[j].myTracerNum=-1 then
        if (pos('_abd',processes[i].input[j].tracer)<=0) and (pos('_mass',processes[i].input[j].tracer)<=0) then //prevent error message being displayed if input is a stage-resolving vector tracer
          s:=s+'Process "'+processes[i].name+'" has input of tracer "'+processes[i].input[j].tracer+'", which is not defined.'+chr(13);

      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
        if lowercase(processes[i].output[j].tracer) = lowercase(tracers[k].name) then
          processes[i].output[j].myTracerNum:=k;
      if processes[i].output[j].myTracerNum=-1 then
        if (pos('_abd',processes[i].output[j].tracer)<=0) and (pos('_mass',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.'+chr(13);
        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+'/'+cellheightTimesDensity;
              if (tracers[k].vertLoc=0) and (processes[ii].vertLoc=2) then
                ratestring:=ratestring+'/'+cellheightTimesDensity;
              if (tracers[k].vertLoc=3) and (processes[ii].vertLoc=0) then
                ratestring:=ratestring+'*'+cellheightTimesDensity;
              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+'/'+cellheightTimesDensity;
              if (tracers[k].vertLoc=0) and (processes[ii].vertLoc=2) then
                ratestring:=ratestring+'/'+cellheightTimesDensity;
              if (tracers[k].vertLoc=3) and (processes[ii].vertLoc=0) then
                ratestring:=ratestring+'*'+cellheightTimesDensity;
              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;
  for i:=length(auxiliaries)-1 downto 0 do
    if ((auxiliaries[i].isZIntegral=1) and (auxiliaries[i].calcAfterProcesses=0)) or (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
           auxiliaries[j].myCalcBeforeZIntegral:=1;
      end;
    end;
  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;
  //finally, 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);


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

end.
