program migratory_flight; {version 2.1 May 1998}

uses Crt, Turbo3, Dos;

const deffile='migrate.def';
const datfile='migrate.dat';
const outfile='migrate.out';
const spffile='migrate.spf';



const g=9.81;                                     {acceleration due to gravity}

const  Windows    = 5;
       Wtab       : array [0..Windows,1..4] of Integer           { X0,Y0,X1,Y1}
                  = ( ( 1,  1,  80, 15),
                      ( 2,  3,  12, 18),
                      (15,  6,  53, 18),
                      (56,  3,  78, 18),
                      ( 2, 21,  78, 23),
                      (15,  3,  53,  3) );
type name=string[4];

type line=string[60];

var varval: array [1..100] of real;
    varnam: array [1..100] of name;
    varlab: array [1..100] of line;
    VFE: array [1..100] of char;
    inp,out,pri: array [1..100] of boolean;
    ifil,ofil:text;
    infil,outfil,eachstep,batch,fullscreen,stepwise,SPSSformat:boolean;
    number:integer;
    accuracy:real;
    cmdline:string;

Procedure Exist (OutIn: word;FileName: line);
var
 Fil: file;
 Present: boolean;
begin
 Assign(Fil, FileName);
 {$I-}
 if OutIn=0 then
  Rewrite(Fil);
 if OutIn=1 then
  Reset(Fil);
 {$I+}
 Present:=(IOresult=0);
 if present then
   close(Fil)
  else
   begin
    writeln;
    writeln(' *** FATAL ERROR ***');
    writeln;
    write(FileName:1);
    if OutIn=1 then
     writeln(' does not exist or is in use by another program');
    if OutIn=0 then
     writeln(' cannot be overwritten. Probably in use by another program');
    writeln;
    writeln(' Program halted');
    halt;
   end;
end;

procedure Frame(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY: Integer);
var  I : Integer;

begin {Frame}
  GotoXY(UpperLeftX, UpperLeftY);
  Write(chr(218));
  for I := (UpperLeftX + 1) to (LowerRightX - 1) do
  begin
    Write(chr(196));
  end;
  Write(chr(191));
  for I := (UpperLeftY + 1) to (LowerRightY - 1) do
  begin
    GotoXY(UpperLeftX , I);  Write(chr(179));
    GotoXY(LowerRightX, I);  Write(chr(179));
  end;
  GotoXY(UpperLeftX, LowerRightY);
  Write(chr(192));
  for I := (UpperLeftX + 1) to (LowerRightX - 1) do
  begin
    Write(chr(196));
  end;
  Write(chr(217));
end; {Frame}


procedure MakeScreen;
var I:integer;
begin
  ClrScr;
 for I := 1 to Windows do
  begin
    Frame(Wtab[I,1] - 1, Wtab[I,2] - 1, Wtab[I,3] + 1, Wtab[I,4] + 1);
    GotoXY(Wtab[I,1],Wtab[I,2]-1);
    Case I of
      1:Write('Flying');
      2:Write('Underway: on outfile');
      3:Write('nr. MFDE MFDW MFRE MFRW');
      4:Write('Flying: final estimates');
      5:Write('Underway: reserves & clock')
     end;
  end;
end; {MakeScreen}


procedure SelectWindow(Win: Integer);
begin
  Window(Wtab[Win,1], Wtab[Win,2], Wtab[Win,3], Wtab[Win,4]);
  case win of
    0: textcolor(lightred);
    1: textcolor(yellow);
    2: textcolor(lightgreen);
    3: textcolor(cyan);
    4: textcolor(lightred);
    5: textcolor(LightGray);
   end;
  textbackground(blue);
end; {SelectWindow}


procedure preset;
var n:integer;
begin
 clrscr;
 number:=0;
 for n:=1 to 100 do
  begin
   varval[n]:=-999;
   VFE[n]:='A';
   pri[n]:=false;
  end;
end;

procedure read_variable_defenitions;
var fil:text;
    n:integer;
    c1,c2,c3,c4,c5,c6,c7,c8:char;

begin
 exist(1,deffile);
 assign(fil,deffile);
 reset(fil);
 while not eof(fil) do
  begin
   repeat                               {TO SKIP COMMENT LINES}
     read(fil,n);
     if n=0 then
      readln(fil);
    until n>0;
   read(fil,c1,varnam[n],c1,c1,c2,c3,c3,c3,c3,c3,c4,VFE[n]);
   if c1='y' then inp[n]:=true else inp[n]:=false;
   if c2='y' then out[n]:=true else out[n]:=false;
   if c3='y' then pri[n]:=true else pri[n]:=false;
   case VFE[n] of
     'v':VFE[n]:='V';
     'f':VFE[n]:='F';
     'e':VFE[n]:='E';
    end;
   if (inp[n]) and (VFE[n]='V') then
     read(fil,varval[n],c1)
    else
     read(fil,c1,c2,c3,c4,c5,c6,c7,c8);
   readln(fil,varlab[n]);
  end;
 close(fil);
end;

procedure variable_overview;

 procedure subset(i,o:boolean);
 var j,n:integer;
     key:char;
  begin
   outfil:=false;
   window(1,3,80,23);
   j:=0;
   for n:=1 to 100 do
    if (VFE[n]<>'A') and (inp[n]=i) and (out[n]=o) then
     begin
      if j=0 then
       begin
        TextColor(LightGreen);
        writeln('V=default value   F=input from file  E=according to equation');
        writeln('var# name VFE   Value out   variable_label');
        TextColor(White);
       end;
      j:=j+1;
      write(n:3,varnam[n]:6,VFE[n]:3);
      if VFE[n]='V' then
        write(varval[n]:10:4)
       else
        write('          ');
      if pri[n] then
        begin
         write(' Y ');
         outfil:=true;
        end
       else
        write(' N ');
      writeln(varlab[n]);
      if j/17=round(j/17) then
       begin
        writeln;
        TextColor(lightred);
        write('press any key to continue.....');
        TextColor(white);
        if not batch then
          read(kbd,key);
        j:=0;
        clrscr;
       end;
     end;
   if j>0 then
    begin
     writeln;
     TextColor(lightred);
     write('press any key to continue.....');
     TextColor(white);
     if not batch then
       read(kbd,key);
    end;
   window(1,1,80,25);
   clrscr;
  end;

begin
 TextColor(yellow);
 writeln('VARIABLES OF WHICH INPUT IS OBLIGATORY');
 TextColor(white);
 subset(true,false);
 TextColor(yellow);
 writeln('VARIABLES OF WHICH INPUT IS FACULTATIVE');
 TextColor(white);
 writeln('when VFE=E the calculated value by the model is further used');
 subset(true,true);
 TextColor(yellow);
 writeln('VARIABLES THAT ARE OUTPUT OF THE MODEL');
 TextColor(white);
 subset(false,true);
end;

procedure open_infile;
var n:integer;
begin
 TextColor(Yellow);
 writeln('VARIABLES READ FROM FILE IN THE FOLLOWING ORDER ARE');
 TextColor(white);
 writeln;
 infil:=false;
 for n:=1 to 100 do
  if VFE[n]='F' then
   begin
    writeln(n:3,varnam[n]:6,'  ',varlab[n]);
    infil:=true;
   end;
 if infil then
   begin
    writeln;
    write('VALUES ARE EXPECTED ON FILE ');
    textcolor(yellow);
    writeln(datfile);
    textcolor(white);
    exist(1,datfile);
    assign(ifil,datfile);
    reset(ifil);
   end
  else
   writeln(' NO VARIABLES ARE READ FROM FILE ! ');
end;

Procedure MakeSpssHeader;
var n,liner,last:integer;
begin
 liner:=1;
 writeln(ofil,'* SPSS JOB HEADER MIGRATORY FLIGHT SIMULATION.');
 writeln(ofil,'DATA LIST FILE=''',outfile,'''/number 1-10');
 for n:=1 to 100 do
  if pri[n] then
   last:=n;
 for n:=1 to last do
   if pri[n] then
    begin
     liner:=liner+1;
     write(ofil,' ',varnam[n]:4,' ',((liner-1)*10+1):1,'-',(liner*10):1);
     if liner=7 then
      begin
       liner:=0;
       if n=last then
         writeln(ofil,'.')
        else
         writeln(ofil,'/');
      end;
    end;
 if liner<>0 then
  writeln(ofil,'.');
 writeln(ofil,'VARIABLE LABELS  number ''program input case number''/');
 for n:=1 to last do
   if pri[n] then
    begin
     write(ofil,'  ',varnam[n]:4,' ''',varlab[n],'''');
     if n=last then
       writeln(ofil,'.')
      else
       writeln(ofil,'/');
    end;
end;
  
Procedure MakeExcelHeader;
var n:integer;
begin
 write(ofil,'     CASE ');
 for n:=1 to 100 do
  if pri[n] then
   write(ofil,varnam[n]:10);
 writeln(ofil);
end;


procedure OfilWriting;                      {output for SPSS/EXCEL to file ofil}
var n,liner,last:integer;
begin
 liner:=1;
 write(ofil,number:10);
 for n:=1 to 100 do
   if pri[n] then
    begin
     liner:=liner+1;
     if round(varval[n])=-999 then
      write(ofil,'          ')
     else
      write(ofil,varval[n]:10:4);
     if liner=7 then
      begin
       liner:=0;
       if SPSSformat then
        writeln(ofil)
      end;
    end;
 if (SPSSformat) and (liner<>0) then
  writeln(ofil);
 if SPSSformat=false then
  writeln(ofil);
end;


procedure open_outfile;
var key:char;
begin
 if outfil then
   begin
    writeln;
    write('SPSS HEADER IS WRITTEN TO ');
    TextColor(yellow);
    write(spffile);
    TextColor(white);
    exist(0,spffile);
    assign(ofil,spffile);
    rewrite(ofil);
    MakeSpssHeader;
    Close(ofil);
    writeln('   READY');
    writeln;
    write('DATA ARE WRITTEN TO FILE ');
    TextColor(yellow);
    writeln(outfile);
    TextColor(white);
    exist(0,outfile);
    assign(ofil,outfile);
    rewrite(ofil);
    writeln;
    write(
    'Should each step be written to the outputfile or the final result only (S/F)');
    if not batch then
      begin
       repeat
        read(kbd,key);
       until key in ['S','s','F','f'];
      end
     else
      key:=cmdline[1];
    writeln;
    if key in ['S','s'] then
      begin
       write(' Each step? OK!  ');
       EachStep:=true;
      end
     else
      begin
       write('    Only finals? ... OK!');
       Eachstep:=false;
      end;
    writeln;
    write(
    'Use output format suitable for importing in SPSS or EXCEL (S/E)');
    if not batch then
      begin
       repeat
        read(kbd,key);
       until key in ['S','s','E','e'];
      end
     else
      key:=cmdline[2];
    writeln;
    if key in ['S','s'] then
      begin
       write('    SPSS format? ... OK!  ');
       SPSSformat:=true;
      end
     else
      begin
       write('    EXCEL format? ... OK!');
       SPSSformat:=false;
       MakeExcelHeader;
      end;
   end
  else
   Eachstep:=false;
 writeln;
 textcolor(lightred);
 write('                          press any key to continue.....');
 textcolor(white);
 if not batch then
   read(kbd,key);
end;


procedure tailwind(V:real); {if MD is entered as a negative value, this procedure}
const pi=3.141592;    {does not correct for drift!, the absolute value of MD}
var alpha,wortel:real;   {is used in the calculations; the negative sign has}
begin                    {thus no value!}
 if VFE[4]='E' then
  begin
   alpha:=(varval[7]+180-abs(varval[8]))*2*pi/360;
   wortel:=sqr(V)-sqr(varval[6]*sin(alpha));
   if wortel<0 then
    wortel:=0-sqrt(0-wortel)
   else
    wortel:=sqrt(wortel);
  if varval[8]<0 then       {if MD negative then no drift correction}
    wortel:=V;
  varval[4]:=varval[6]*cos(alpha)+wortel-V;
 end;
end;



procedure pennycuick(mass:real);
var Sb,Cdb,A,Sd,Pam,Ppro,Pmet,x2,vmppril,Pmin,Dpar,Ppar,rf,Pind,V,p,Vextra:real;
    leanmass: real;
    I,Vlow,Vupp:integer;
    key:char;
begin
                                                            {body frontal area}
 Sb:=(8.13*0.001)*exp(0.666*ln(mass));
                                                         {body Reynolds number}
 if VFE[27]='E' then
   begin
    Cdb:=125000.0*sqrt(mass);
    if Cdb<50000.0 then
      Cdb:=0.4
     else
      if Cdb>200000.0 then
        Cdb:=0.25
       else
        Cdb:=1.57-(0.108*ln(Cdb));
   end
  else
   Cdb:=varval[27];

                                                   {Equivalent flat plate area}
 A:=Sb*Cdb;
                                                                    {disc area}
 Sd:=3.142*sqr(varval[14])/4;
                                                       {absolute minimum power}
 Pam:=0.877*exp(0.75*ln(varval[22]))*exp(1.5*ln(mass*g))*exp(0.25*ln(A))
     /( sqrt(varval[20])*exp(0.75*ln(Sd)) );
                                                                {profile power}
 Ppro:=varval[23]*Pam;
                        {Find BMR and combine with power conversion efficiency}
                        {BMR calculated over lean body mass}
 leanmass:=varval[11]*(1-varval[10]);
 if round(varval[15])=1 then
   Pmet:=varval[21]*6.25*exp(0.724*ln(leanmass))
  else
   Pmet:=varval[21]*3.79*exp(0.723*ln(leanmass));
 x2:=Pmet/Pam;
                                                  {preliminary estimate of vmp}
 vmppril:=(0.76* exp(0.25*ln(varval[22])) *sqrt(mass*g))/
          (sqrt(varval[20])* exp(0.25*ln(Sd*A)) );
                                                                  {lower limit}
 Vlow:=10*trunc(0.8*vmppril);
                                                                  {upper limit}
 if (varval[23]+x2)<2.5 then
   Vupp:=10*(trunc(1.9*vmppril)+1)
  else
   Vupp:=10*(trunc(2.3*vmppril)+1);
 if VFE[32]='E' then               {check if user defined value for Vmr exists}
   varval[32]:=0
  else
   begin
    Vlow:=round(varval[32]*10);
    Vupp:=round(varval[32]*10);
   end;
 varval[33]:=0;
 varval[34]:=0.00000000001;
 Pmin:=0;                          if fullscreen then
                                    begin
                                     SelectWindow(1);
                                     Write('Va,m/s  P,W');
                                    end;
 for I:=Vlow to Vupp do
  begin
   V:=I/10;
                                                                {parasite drag}
   Dpar:=varval[20]*sqr(V)*A/2;
                                                               {parasite power}
   Ppar:=Dpar*V;
                                            {resultant of weight and body drag}
   rf:=sqrt(sqr(mass*g)+sqr(Dpar));
                                                                {induced power}
   Pind:=varval[22]*sqr(rf)/(2*Sd*V*varval[20]);
   p:=varval[24]*(Ppar+Pind+Ppro+Pmet);
   if (I=Vlow) or (p<Pmin) then
    begin
     varval[33]:=V;
     Pmin:=p;
    end;
   if round(varval[26])=1 then
      Vextra:=0
     else
      begin
       tailwind(V);
       Vextra:=varval[4];
      end;
   if (I=Vlow) or ((V+Vextra)/p>(varval[32]+Vextra)/varval[34]) then
    begin
     varval[32]:=V;
     varval[34]:=p;
    end;
                            if (trunc(V)=V) and (fullscreen) then
                              writeln(V:2:0,p:7:2);
  end;
                            if fullscreen then
                             begin
                              writeln('         ');
                              write('         ');
                                 {Extra check possibility on Pennyquicks model}
                              SelectWindow(4);
                              write('mass',mass:7:4);
                              write(' kg      b  ',varval[14]:8:3);
                              write(' m       Sd ',Sd:8:3);
                              writeln(' m2      Sb ',Sb:8:5,' m2 ');
                              write('Cdb',Cdb:8:3);
                              write('         A  ',A:8:5);
                              write(' m2      Pmet',Pmet:7:3);
                              writeln(' W       X2 ',x2:8:3,'  ');
                              write('Vmp',varval[33]:8:3);
                              write(' m/s     Vmr',varval[32]:8:3);
                              write(' m/s     Pmin',Pmin:7:3);
                              write(' W       Pmr',varval[34]:8:3,' W ');
                             end;
end;


function wvd(T,RH,press:real):real;               {water vapor density g/l STPD}
   {according to F.W. Murray 1967   [in g/L]}
   {On the computation of saturation vapor pressure}
   {J. Applied Meteorology 6: 203-204}
 var wvdunstd: real;
     correction: real;
     alpha: real;
 begin
   alpha:=22.414/18;
   correction:=1013*(273.15+T)/(Press*273.15);
   wvdunstd:=RH/100*exp((16.78*T-116.9)/(T+237.3))*18/(8.31432*(273.2+T));
   wvd:=wvdunstd*correction/(1-alpha*wvdunstd*correction);
 end;

function AirDensity(P,T,RH: real):real;                   {Air density kg/m^3}
 var E,Dw,Da:real;
begin
  {evaporation pressure in hPa}
     E:=6.11*exp(7.5*T/(237.7+T)*ln(10))*RH/100;
  {density of water vapour}
     Dw:=P*100/((T+273.15)*461.5);
  {density of dry air}
     Da:=P*100/((T+273.15)*287);
  AirDensity:=E/P*Dw+(1-E/P)*Da;
end;

procedure model;                                 {here the real stuff happens}
var maxenergy,maxwater:boolean;
    mass,MFDE,MFDW,MFRE,MFRW,oldvv75,oldvv76:real;
    time,n:integer;
    first:boolean;
    key:char;

                    procedure topline;   {procedure to actualize command line}
                    begin
                     SelectWindow(0);
                     write(' Command keys:  Quit');
                     if stepwise then
                       write('  Continue')
                      else
                       write('  Stepwise');
                     if fullscreen then
                       write('  Partial screen output')
                      else
                       write('  Full screen output   ');
                     if stepwise then
                       begin
                        textcolor(blink);
                        write('      Press space...');
                       end
                      else
                        write('                    ');
                    end;


                    procedure MustIStopNow;   {to invoke immediate and total}
                     var Ch:char;               {break on the user's request}
                    begin                    
                     if keypressed then
                      begin
                       Ch:=readkey;
                       if Ch in ['Q','q'] then
                        begin
                         stepwise:=false;
                         maxenergy:=true;
                         maxwater:=true;
                         if infil then
                           while not eof(ifil) do
                            readln(ifil);
                        end;
                       if Ch in ['S','s'] then
                        stepwise:=true;
                       if Ch in ['C','c'] then
                        stepwise:=false;
                       if Ch in ['F','f'] then
                        fullscreen:=true;
                       if Ch in ['P','p'] then
                        fullscreen:=false;
                       topline;
                      end;
                    end;

                                             {start of the model procedure}
begin                                           {initializing a few things}
 number:=number+1;
 if number=1 then
  begin
   if batch then
    begin
     fullscreen:=false;
     stepwise:=false;
    end
   else
    begin
     fullscreen:=true;
     stepwise:=true;
    end;
   topline;
  end;
 maxwater:=false;
 maxenergy:=false;
 accuracy:=round(varval[77])*60;
 time:=0;
 for n:=71 to 74 do                           {set MFDE MFDW MFRE MFRW to zero}
  varval[n]:=0;
 varval[75]:=varval[11]*varval[10]*varval[41]*1000.0;  {kJ}
 varval[76]:=varval[11]*varval[12]*varval[13]*1000;        {g}
   {set a few parameters that are not affected by body mass loss during flight} 
                                                   {air pressure, mbar=hpascal}
 if VFE[5]='E' then                        {check if user defined value exists}
   varval[5]:=(1013*exp(4.256*ln(1-0.00000688*varval[1]/0.3048)));
                                                                  {air density}
 if VFE[20]='E' then                       {check if user defined value exists}
   varval[20]:=AirDensity(varval[5],varval[2],varval[3]);
                                                    {exhaled air temperature C}
 if VFE[58]='E' then                       {check if user defined value exists}
   varval[58]:=23.0+0.43*varval[2];
 mass:=varval[11];
                              {we start flying and flying until we hit a limit}
 repeat
  MustIStopNow;                                        {? reacts on keypressed}
 
              {POWER FOR FLIGHT AND METABOLIC WATER PRODUCTION}
                                              {calculate flight speed and power}
  pennycuick(mass);
                                                       {maximum lift drag ratio}
  varval[31]:=mass*g*varval[32]/varval[34];
                                        {aerobic or chemical power input, Watt}
  if VFE[35]='E' then                      {check if user defined value exists}
    varval[35]:=varval[34]/varval[21];
                                                        {fuel consumption mg/h}
  varval[42]:=varval[35]*3600/varval[41];
                                              {metabolic water production mg/h}
  varval[44]:=varval[42]*varval[43];

                   {MINIMUM EVAPORATIVE WATER LOSS MODEL}
                                          {rate of oxygen consumption STP ml/h}
  varval[52]:=varval[35]*3600/varval[51];
                                            {pulmonary ventilation volume ml/h}
  varval[57]:=varval[52]/(0.2095*varval[54]);
                                      {respiratory evaporative water loss mg/h}
  varval[59]:=(varval[57]+varval[52]*(varval[53]-1))
                   *wvd(varval[58],100,varval[5])
                   -varval[57]*wvd(varval[2],varval[3],varval[5]);
                                        {cutaneous evaporative water loss mg/h}
  if VFE[60]='E' then                      {check if user defined value exists}
                                 {calculate CWL according to Kvist et al. 1998}
   varval[60]:=1/(1+varval[59]/ 
                      (3.48*3600*exp(0.667*ln(mass*1000))*0.001) );
                                            {total evaporative water loss mg/h}
  if VFE[61]='E' then                      {check if user defined value exists}
   varval[61]:=varval[59]/(1-varval[60]);

                      {EVAPORATIVE HEAT REGULATION MODEL}
                                                     {standard dry conductance} 
  if VFE[56]='E' then                      {check if user defined value exists}
   varval[56]:=0.0112*exp(0.517*ln(mass*1000));
                                            {total evaporative water loss mg/h}
  if VFE[62]='E' then                      {check if user defined value exists}
   begin
    varval[62]:=(  varval[34]*(1/varval[21]-1)-
                   varval[56]*(1+0.26*sqrt(varval[32]))*(varval[19]-varval[9]))
                   /varval[50]*1000*3600;
   end;
                               {NET WATER LOSS}
                                                  {net rate of water loss mg/h}
  if VFE[63]='E' then                      {check if user defined value exists}
   if varval[61]>varval[62] then
     varval[63]:=varval[61]-varval[44]
    else
     varval[63]:=varval[62]-varval[44];
                                             {increase variables for this step}
  tailwind(varval[32]);
  if maxenergy=false then                                     {set MFDE & MFRE}
   if varval[35]<=0 then       {if no change in energy reserves or increasing!}
     begin
      varval[71]:=-999;
      varval[73]:=-999;
     end
    else                             {if appropriate change in energy reserves}
     begin
      varval[71]:=varval[71]+accuracy/3600;
      varval[73]:=varval[73]+(varval[32]+varval[4])*accuracy/1000;
     end;
  if maxwater=false then
     begin
      varval[72]:=varval[72]+accuracy/3600;
      varval[74]:=varval[74]+(varval[32]+varval[4])*accuracy/1000;
     end;
  oldvv75:=varval[75];
  oldvv76:=varval[76];
  varval[75]:=varval[75]-(varval[35]*accuracy)/1000;
  varval[76]:=varval[76]-(varval[63]/1000)*(accuracy/3600);
  if round(varval[18])=1 then
    if (varval[63]<0) and (varval[76]>=varval[11]*varval[12]*varval[13]*1000) then
      mass:=mass-((varval[42]+     0.0  )*(accuracy/3600))/1000000.0
     else
      mass:=mass-((varval[42]+varval[63])*(accuracy/3600))/1000000.0;
   {else: no correction for on flight water and fat loss}
  if varval[76]>varval[11]*varval[12]*varval[13]*1000 then
                     {correct for weight gain due to positive water balance}
      varval[76]:=varval[11]*varval[12]*varval[13]*1000;
  time:=time+1;
  if (varval[75]<0) or (varval[35]<=0) or (mass<=varval[11]*(1-varval[17])) then
    maxenergy:=true;
  if (varval[76]<0) or (mass<=varval[11]*(1-varval[17])) then
    maxwater:=true;
  if ((maxwater) or (maxenergy)) and (varval[79]=2) then
    begin
     maxenergy:=true;
     maxwater:=true;
    end;

  {interpolations for accurate determination of flight ranges and durations}
  if (maxenergy) and (varval[75]<0) and (oldvv75>0) then
   begin
    varval[71]:=varval[71]+varval[75]*accuracy/(oldvv75-varval[75])/3600;
    varval[73]:=varval[73]+varval[75]*accuracy/(oldvv75-varval[75])*(varval[32]+varval[4])/1000;
    if varval[79]=2 then
     begin
      varval[72]:=varval[71];
      varval[74]:=varval[73];
     end;
   end;
  if (maxwater) and (varval[76]<0) and (oldvv76>0) then
   begin
    varval[72]:=varval[72]+varval[76]*accuracy/(oldvv76-varval[76])/3600;
    varval[74]:=varval[74]+varval[76]*accuracy/(oldvv76-varval[76])*(varval[32]+varval[4])/1000;
    if varval[79]=2 then
     begin
      if varval[72]<varval[71] then
        begin
         varval[71]:=varval[72];
         varval[73]:=varval[74];
        end
       else
        begin
         varval[72]:=varval[71];
         varval[74]:=varval[73];
        end
     end;
   end;
                           if fullscreen then
                             begin
                              selectwindow(5);
                              write('Eres ',varval[75]:1:0,'kJ');
                              write('  Wres ',varval[76]:3:1,'g');
                              write('  Time ',time*(accuracy/3600):3:1,'h');
                              Clreol;
                              selectwindow(2);
                              first:=true;
                              for n:=1 to 100 do
                               if pri[n] then
                                begin
                                 if (first) then
                                   write(varnam[n]:4,' ')
                                  else
                                   write('         ',varnam[n]:4,' ');
                                 if round(varval[n])<>-999 then
                                   write(varval[n]:9:4)
                                  else
                                   write('         ');
                                 if first then
                                   first:=false
                                  else
                                   begin
                                    writeln;
                                    first:=true;
                                   end;
                                end;
                             end;
  if outfil then
   begin
    if EachStep then
      OfilWriting;
    if (Eachstep=false) and (((maxenergy) and (maxwater))
                                          or (varval[78]<=time)) then
      OfilWriting;
   end;
  if stepwise then
    repeat until keypressed;
 until (((maxenergy) and (maxwater)) or (varval[78]<=time));
                                        selectwindow(3);
                                        GotoXY(1,Wtab[3,4]-Wtab[3,2]+1);
                                        InsLine;
                                        write(number:2,' ');
                                        for n:=71 to 74 do
                                         if round(varval[n])>=0 then
                                           write(varval[n]:4:0,' ')
                                          else
                                           write('  -- ');
end;

procedure repeated_model;
var n:integer;
begin
 while not eof(ifil) do
  begin
   for n:=1 to 100 do
    if VFE[n]='F' then
      read(ifil,varval[n]);
   readln(ifil);
   model;
  end;
end;

procedure close_files;
begin
 if infil then
  close(ifil);
 close(ofil);
end;



{MAIN PROGRAM MAIN PROGRAM MAIN PROGRAM MAIN PROGRAM MAIN PROGRAM MAIN PROGRAM}

begin
 textbackground(blue);
 textcolor(white);
 preset;
 cmdline:=paramstr(1);

{ cmdline:='f'; }

 if cmdline='' then
   batch:=false
  else
   batch:=true;
 read_variable_defenitions;
 if not batch then
   variable_overview;
 open_infile;
 open_outfile;
 makescreen;
 if infil then
   repeated_model
  else
   model;
 close_files;
 Window(1,1,80,25);
 GotoXY(1,24);
 textbackground(black);
 textcolor(white);
end.