unit erg_pseudo3d;

//special treatment of tracers and processes with vertLoc=3
//defines auxiliary variables and modifies processes

interface

uses classes;

procedure TreatPseudo3d(var s: TStringList);

implementation

uses erg_base, erg_types, sysUtils;

procedure TreatPseudo3d(var s: TStringList);
var
  i, j, k, l, jmax, lmax: Integer;
  found, found2: Boolean;
begin
  //first, seek for tracers with vertLoc=3 to create a zIntegral of their "verticalDistribution" auxiliary variable
  for i:=0 to length(tracers)-1 do
    if tracers[i].vertLoc=3 then
    begin
      found:=false;
      if (tracers[i].verticalDistribution='1') or (tracers[i].verticalDistribution='1.0') then
      begin
        found:=true;
        //we want to create an auxiliary variable unity_zIntegral, unless it already exists
        found2:=false;
        for k:=0 to length(auxiliaries)-1 do
          if lowercase(auxiliaries[k].name) = 'unity_zintegral' then found2:=true;
        if found2=false then
        begin
          setLength(auxiliaries,length(auxiliaries)+1);
          InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
          auxiliaries[length(auxiliaries)-1].name:='unity_zIntegral';
          auxiliaries[length(auxiliaries)-1].formula:='1.0';
          auxiliaries[length(auxiliaries)-1].vertLoc:=1;
          auxiliaries[length(auxiliaries)-1].isZIntegral:=1;
          auxiliaries[length(auxiliaries)-1].description:='vertical integral of density [kg/m2]';
        end;
      end
      else
      begin
        jmax:=length(auxiliaries)-1;
        for j:=0 to jmax do
          if trim(lowercase(auxiliaries[j].name))=trim(lowercase(tracers[i].verticalDistribution)) then
          begin
            found:=true;
            if auxiliaries[j].vertLoc<>0 then
              s.Add('Auxiliary variable '+auxiliaries[j].name+' was chosen for verticalDistribution of tracer '+tracers[i].name+', but has vertLoc/=0')
            else
            begin
              //we want to create a zIntegral of this auxiliary variable, unless it already exists
              found2:=false;
              for k:=0 to length(auxiliaries)-1 do
                if lowercase(auxiliaries[k].name) = lowercase(auxiliaries[j].name)+'_zintegral' then found2:=true;
              if found2=false then
              begin
                setLength(auxiliaries,length(auxiliaries)+1);
                InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
                auxiliaries[length(auxiliaries)-1].name:=auxiliaries[j].name+'_zIntegral';
                auxiliaries[length(auxiliaries)-1].formula:=auxiliaries[j].name;
                auxiliaries[length(auxiliaries)-1].vertLoc:=1;
                auxiliaries[length(auxiliaries)-1].isZIntegral:=1;
                auxiliaries[length(auxiliaries)-1].description:='vertical integral of ('+auxiliaries[j].description+')*density [kg/m2]';
              end;
            end;
          end;
      end;
      if found=false then s.Add('Tracer '+tracers[i].name+' has verticalDistribution='+tracers[i].verticalDistribution+', but this is neither 1.0 nor an auxiliary variable.');
    end;

  //second, seek for processes with vertLoc=3 which have a tracer with vertLoc=3 in their output. They may only have one input tracer.
  for i:=0 to length(processes)-1 do
    if processes[i].vertLoc=3 then
    begin
      //first, seek if a tracer with vertLoc=3 exists in the output
      found:=false;
      for k:=0 to length(processes[i].output)-1 do
        if tracers[processes[i].output[k].myTracerNum].vertLoc=3 then
        begin
          if found then
            s.Add('Process '+processes[i].name+' has output of more than one tracer with vertLoc=3, which is forbidden.');
          found:=true;
          j:=processes[i].output[k].myTracerNum;
        end;
      if found then //Feeding process
      begin
          //this is a feeding process, it may only have one input tracer. Otherwise we get in trouble with stochiometry when the positive-definite scheme limits the process in a certain depth.
          if length(processes[i].input)>1 then
            s.Add('Feeding process '+processes[i].name+' has input of more than one tracer, which is forbidden.')
          else if length(processes[i].input)=1 then
          begin
            k:=processes[i].input[0].myTracerNum;
            //seek for a tracer with vertLoc=1 in input or output. If this exists, process takes place in the bottom cell only.
            found:=false;
            for l:=0 to length(processes[i].output)-1 do
              if tracers[processes[i].output[l].myTracerNum].vertLoc=1 then
                found:=true;
            if (tracers[k].vertLoc=0) and (not found) then
            begin
              //input-tracer has vertLoc=0, no output tracer with vertLoc=1 -> process gets vertLoc=0
              //step 1/4: create auxiliary variable tra1_seen_by_tra2_3d
              lmax:=length(auxiliaries)-1;
              found:=false;
              for l:=0 to lmax do
                if lowercase(auxiliaries[l].name)=lowercase(tracers[k].name)+'_seen_by_'+lowercase(tracers[j].name)+'_3d' then
                  found:=true;
              if found=false then
              begin
                setLength(auxiliaries,length(auxiliaries)+1);
                InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
                auxiliaries[length(auxiliaries)-1].name:=tracers[k].name+'_seen_by_'+tracers[j].name+'_3d';
                auxiliaries[length(auxiliaries)-1].formula:=tracers[k].name+'*'+processes[i].feedingEfficiency;
                auxiliaries[length(auxiliaries)-1].description:=tracers[k].description+' which can be preyed upon by '+tracers[j].description+' [mol/kg]';
              end;
              //step 2/4: create auxiliary variable tra1_seen_by_tra2
              lmax:=length(auxiliaries)-1;
              found:=false;
              for l:=0 to lmax do
                if lowercase(auxiliaries[l].name)=lowercase(tracers[k].name)+'_seen_by_'+lowercase(tracers[j].name) then
                  found:=true;
              if found=false then
              begin
                setLength(auxiliaries,length(auxiliaries)+1);
                InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
                auxiliaries[length(auxiliaries)-1].name:=tracers[k].name+'_seen_by_'+tracers[j].name;
                auxiliaries[length(auxiliaries)-1].formula:=tracers[k].name+'_seen_by_'+tracers[j].name+'_3d';
                auxiliaries[length(auxiliaries)-1].description:='vertical integral of '+tracers[k].description+' which can be preyed upon by '+tracers[j].description+' [mol/m2]';
                auxiliaries[length(auxiliaries)-1].vertLoc:=1;
                auxiliaries[length(auxiliaries)-1].isZIntegral:=1;
              end;
              //step 3/4: set process.vertLoc to 0
              processes[i].vertLoc:=0;
              //step 4/4: convert process rate from mol/m2/day to mol/kg/day
              processes[i].turnover:='('+processes[i].turnover+')*'+tracers[k].name+'_seen_by_'+tracers[j].name+'_3d/max(1.0e-30,'+tracers[k].name+'_seen_by_'+tracers[j].name+')';
            end
            else if (tracers[k].vertLoc=1) or found then
            begin
              //input-tracer or one of the output tracers has vertLoc=1
              //step 1/2: create auxiliary variable tra1_seen_by_tra2
              lmax:=length(auxiliaries)-1;
              found:=false;
              for l:=0 to lmax do
                if lowercase(auxiliaries[l].name)=lowercase(tracers[k].name)+'_seen_by_'+lowercase(tracers[j].name) then
                  found:=true;
              if found=false then
              begin
                setLength(auxiliaries,length(auxiliaries)+1);
                InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
                auxiliaries[length(auxiliaries)-1].name:=tracers[k].name+'_seen_by_'+tracers[j].name;
                auxiliaries[length(auxiliaries)-1].formula:=tracers[k].name+'*'+processes[i].feedingEfficiency;
                auxiliaries[length(auxiliaries)-1].description:=tracers[k].description+' which can be preyed upon by '+tracers[j].description+' [mol/m2]';
                auxiliaries[length(auxiliaries)-1].vertLoc:=1;
              end;
              //step 2/2: set process.vertLoc to 1
              processes[i].vertLoc:=1;
            end
            else if tracers[k].vertLoc=3 then
            begin
              //input-tracer has vertLoc=3
              if length(processes[i].output)>1 then
                s.Add('In process '+processes[i].name+', an vertLoc=3 tracer exists both as input and output. Therefore, only one input and output tracer are allowed.')
              else
              begin
                //step 1/4: create auxiliary variable tra1_seen_by_tra2_3d
                lmax:=length(auxiliaries)-1;
                found:=false;
                for l:=0 to lmax do
                  if lowercase(auxiliaries[l].name)=lowercase(tracers[k].name)+'_seen_by_'+lowercase(tracers[j].name)+'_3d' then
                    found:=true;
                if found=false then
                begin
                  setLength(auxiliaries,length(auxiliaries)+1);
                  InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
                  auxiliaries[length(auxiliaries)-1].name:=tracers[k].name+'_seen_by_'+tracers[j].name+'_3d';
                  auxiliaries[length(auxiliaries)-1].formula:=tracers[k].verticalDistribution+'*'+processes[i].feedingEfficiency;
                  auxiliaries[length(auxiliaries)-1].description:='vertical distribution of '+tracers[k].description+' which can be preyed upon by '+tracers[j].description+' [1]';
                end;
                //step 2/4: create auxiliary variable tra1_seen_by_tra2_3d_zintegral
                lmax:=length(auxiliaries)-1;
                found:=false;
                for l:=0 to lmax do
                  if lowercase(auxiliaries[l].name)=lowercase(tracers[k].name)+'_seen_by_'+lowercase(tracers[j].name)+'_3d_zintegral' then
                    found:=true;
                if found=false then
                begin
                  setLength(auxiliaries,length(auxiliaries)+1);
                  InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
                  auxiliaries[length(auxiliaries)-1].name:=tracers[k].name+'_seen_by_'+tracers[j].name+'_3d_zIntegral';
                  auxiliaries[length(auxiliaries)-1].formula:=tracers[k].name+'_seen_by_'+tracers[j].name+'_3d';
                  auxiliaries[length(auxiliaries)-1].description:='vertical integral of vertical distribution of '+tracers[k].description+' which can be preyed upon by '+tracers[j].description+', times density [kg/m2]';
                  auxiliaries[length(auxiliaries)-1].vertLoc:=1;
                  auxiliaries[length(auxiliaries)-1].isZIntegral:=1;
                end;
                //step 3/4: create auxiliary variable tra1_seen_by_tra2
                lmax:=length(auxiliaries)-1;
                found:=false;
                for l:=0 to lmax do
                  if lowercase(auxiliaries[l].name)=lowercase(tracers[k].name)+'_seen_by_'+lowercase(tracers[j].name) then
                    found:=true;
                if found=false then
                begin
                  setLength(auxiliaries,length(auxiliaries)+1);
                  InitErgAuxiliary(auxiliaries[length(auxiliaries)-1]);
                  auxiliaries[length(auxiliaries)-1].name:=tracers[k].name+'_seen_by_'+tracers[j].name;
                  if (tracers[k].verticalDistribution='1') or (tracers[k].verticalDistribution='1.0') then
                    auxiliaries[length(auxiliaries)-1].formula:=tracers[k].name+'*'+tracers[k].name+'_seen_by_'+tracers[j].name+'_3d_zIntegral/max(unity_zIntegral,1.0e-30)'
                  else
                    auxiliaries[length(auxiliaries)-1].formula:=tracers[k].name+'*'+tracers[k].name+'_seen_by_'+tracers[j].name+'_3d_zIntegral/max('+tracers[k].verticalDistribution+'_zIntegral,1.0e-30)';
                  auxiliaries[length(auxiliaries)-1].description:=tracers[k].description+' which can be preyed upon by '+tracers[j].description+' [mol/m2]';
                  auxiliaries[length(auxiliaries)-1].vertLoc:=1;
                end;
                //step 4/4: set process.vertLoc to 1
                processes[i].vertLoc:=1;
              end;
            end;
          end
          else
          begin
            //no input tracer at all
            processes[i].vertLoc:=1;
          end;
      end;
    end;






  //third, seek for processes with vertLoc=3 which have a tracer with vertLoc=3 in their input, but not in their output
  for i:=0 to length(processes)-1 do
    if processes[i].vertLoc=3 then
    begin
      //seek vertLoc=3 tracer in the input
      found:=false;
      for j:=0 to length(processes[i].input)-1 do
      begin
        if tracers[processes[i].input[j].myTracerNum].vertLoc=3 then
        begin
          if found=false then found:=true
          else
            s.Add('Process '+processes[i].name+' contains more than one input tracer with vertLoc=3, which is prohibited');
          k:=processes[i].input[j].myTracerNum;
        end;
      end;
      if found then //an input tracer with vertLoc=3 was found, it is tracer k
      begin
        //check if another tracer with vertLoc=3 exists in the output - this must not be the case
        found:=false;
        for j:=0 to length(processes[i].output)-1 do
          if tracers[processes[i].output[j].myTracerNum].vertLoc=3 then found:=true;
        if found then
          s.Add('Process '+processes[i].name+' contains an output tracer with vertLoc=3 and also other output tracers. This is prohibited.')
        else
        begin
          //now check if a tracer with vertLoc=1 occurs in the input or in the output
          found:=false;
          for j:=0 to length(processes[i].input)-1 do
            if tracers[processes[i].input[j].myTracerNum].vertLoc=1 then found:=true;
          for j:=0 to length(processes[i].output)-1 do
            if tracers[processes[i].output[j].myTracerNum].vertLoc=1 then found:=true;
          if found then
          begin
            //this process will only happen at the bottom
            processes[i].vertLoc:=1;
          end
          else
          begin
            //maybe no output tracer exists at all? Then process will only happen at the bottom
            if length(processes[i].output)=0 then
              processes[i].vertLoc:=1
            else
            begin
              //this process has output of tracers with vertLoc=0 only and thus will happen in the whole water column
              processes[i].vertLoc:=0;
              if (tracers[k].verticalDistribution='1') or (tracers[k].verticalDistribution='1.0') then
                processes[i].turnover:='('+processes[i].turnover+')*'+tracers[k].verticalDistribution+'/unity_zIntegral'
              else
                processes[i].turnover:='('+processes[i].turnover+')*'+tracers[k].verticalDistribution+'/'+tracers[k].verticalDistribution+'_zIntegral';
            end;
          end;
        end;
      end;
    end;
end;

end.
