program Burstsimulator; // Written by Hannes Tammet
{$APPTYPE CONSOLE}
{$R+}                 
uses SysUtils;

const original = 'HT20060912';
       revised = 'HT20110315';
      revision = original + ' revised ' + revised;
      lastyear = 2011;
     lastmonth = 6;
       lastday = 30;
{The program will be usable until the defined last date.
 The definition lastyear = 9999 cancels the restriction.}

{The program and its parts can be freely used and modified.
 Please give reference to the papers by Tammet and Kulmala:
 1. Simulation tool for atmospheric aerosol nucleation bursts,
    J. Aerosol Sci., 2006, vol. 36, pp. 173-196.
 2. Simulating of aerosol nucleation bursts in conifer forest,
    Boreal Environment Research, 2006, vol. 12, pp. 421-430.

 When modifying the program, at first a revision comment must
 be added and the revision constant must be corrected above.

 Two temporary corrections will facilitate the debugging:
  1) change $R- to $R+ in the program header,
  2) delete "//" before "try" and "except end" 
     (use search to find these words)}

function Mobility                         {  air     nitrogen }
 {velocity/force} (GasMass         {amu}, { 28.96    28.02    }
 {  (m/s) / fN  }  Polarizability  {nm3}, {  0.00171  0.00174 }
                   VisCon1         {nm},  {  0.3036   0.2996  }
 {JAS26, 1995}     VisCon2         {K},   { 44       40       }
 {pp. 459-475}     VisCon3,               {  0.8      0.7     }
 {C: H.Tammet}     Pressure        {mb},
                   Temperature     {K},
                   ParticleDensity {g cm-3, for cluster ions typically 2.08},
                   ParticleCharge  {e, for cluster ions 1},
                   MassDiameter    {nm} : real) : real;
 {NB: Multiply with 1.602 to get electrical mobility of
      a single charged cluster or particle in cm2V-1s-1}
                    
  function Omega11 (x : real) : real; {*(1,1)*(T*) for (*-4) potential}
    var p, q : real;                    {and elastic-specular collisions}
    begin
      if x > 1 then Omega11 := 1 + 0.106 / x + 0.263 / exp ((4/3) * ln (x))
      else begin p := sqrt (x); q := sqrt (p);
        Omega11 := 1.4691 / p - 0.341 / q + 0.181 * x * q + 0.059 end;
    end;

  const a = 1.2; b = 0.5; c = 1; {the slip factor coefficients}
        ExtraDistance = 0.115 {nm}; TransitionDiameter = 2.48 {nm};
  var   GasDiameter, MeanVelocity, Viscosity, FreePath, DipolEffect,
        DeltaTemperature, CheckMark, ParticleMass, CollisionDistance,
        Kn, Omega, s, x, y : real;

  begin
    Viscosity {microPa s} := 0.02713 * sqrt (GasMass * Temperature) /
      sqr (VisCon1 * (1 + exp (VisCon3 * ln (VisCon2 / Temperature))));
    MeanVelocity {m/s} := 145.5 * sqrt (Temperature / GasMass);
    FreePath  {nm} := (166251 * Viscosity * Temperature) /
                      (GasMass * Pressure * MeanVelocity);
    ParticleMass {amu} := 315.3 * ParticleDensity *
                                  exp (3 * ln (MassDiameter));
    DeltaTemperature := Temperature;
    repeat
      CheckMark := DeltaTemperature;
      GasDiameter {nm} := VisCon1 *
        (1 + exp (VisCon3 * ln (VisCon2 / DeltaTemperature)));
      CollisionDistance {nm} := MassDiameter / 2 + ExtraDistance +
                                 GasDiameter / 2;
      DipolEffect := 8355 * sqr (ParticleCharge) * Polarizability /
                            sqr (sqr (CollisionDistance));
                           {NB: erratum (ParticleCharge) in JAS is corrected}
      DeltaTemperature := Temperature + DipolEffect;
    until abs (CheckMark - DeltaTemperature) < 0.01;
    if ParticleCharge = 0 then Omega := 1
                          else Omega := Omega11 (Temperature / DipolEffect);
    Kn := FreePath / CollisionDistance;
    if Kn < 0.03 {underflow safe} then y := 0 else y := exp (- c / Kn);
                 {NB! erratum in JAS (y := 1) corrected!}
    x := (273.15 / DeltaTemperature) *
         exp (3 * ln (TransitionDiameter / MassDiameter));
    if x > 30 {overflow safe} then s := 1
      else if x > 0.001
      then s := 1 + exp (x) * sqr (x / (exp (x) - 1)) * (2.25 / (a + b) - 1)
      else {underflow safe}  s := 1 + (2.25 / (a + b) - 1);
    Mobility := ((2.25 / (a + b)) / (Omega + s - 1))  *
                sqrt (1 + GasMass / ParticleMass) *
                (1 + Kn * (a + b * y)) /
                (6 * PI * Viscosity * CollisionDistance);
  end;

function MassDiameter             {Air environment}
{Uses external function "Mobility"}                                  
         {nm}    (Pressure        {mb},
                  Temperature     {K},
                  ParticleDensity {g cm-3, for cluster ions typically 2.08},
                  ParticleCharge  {e, for cluster ions 1},
                  MechMobility    {m fN-1 s-1} : real) : real;
                 {MechMobility = 0.624 * Z (cm2V-1s-1) / q (e)}
  var c, d, test : real;
               n : integer;
  begin
    c := 300; n := 0;
    repeat
      n := n + 1;
      d := (0.6 + sqrt (0.36 + 200 * c * MechMobility)) /
           (c * MechMobility) - 0.3; 
      test := Mobility (28.96, 0.00171, 0.3036, 44, 0.8,  
                        Pressure, Temperature,
                        ParticleDensity, ParticleCharge, d);
      c := (1.2 / (d + 0.3) + 200 / sqr (d + 0.3)) / test;
    until (abs (test / MechMobility - 1) < 0.0001) or (n = 99);
    if n < 99 then MassDiameter := d else MassDiameter := 0;
  end;

{Ion-particle attachment or coagulation coefficient}
function beta (q, {number of charges on particle, attracting -, repelling +}
{cm3/s}       di, {diameter of ion : nm}
            gcm3, {ion density : g cm-3, typically 2.08}
              dp, {diameter of particle : nm}
               T, {temperature : K}
               p  {pressure : mb} : real) : real;
   {Uses external function "Mobility" (B: 1e15 m s-1 N-1),
    diameters of ions of mobility 1.36 and 1.56 cm2/Vs,
       are in standard conditions 0.79 and 0.70 nm}
   var d, pi2Di, pi2Dp, x, f, y : real;
   begin
      pi2Di := 8.674e-8 * T * mobility (28.96, 0.00171, 0.3036, 44, 0.8,
       {2*pi*D : SI, D = kTB}           p, T, gcm3, 1, di);
      pi2Dp := 8.674e-8 * T * mobility (28.96, 0.00171, 0.3036, 44, 0.8,
       {2*pi*D : SI, D = kTB}           p, T, gcm3, abs (q), dp);
      d := dp + di + 0.23; {nm}
      x := 2 * q * (1.671e4 / T) / d;
      y := 1 - 1 / (1 + 0.3 * q * (q - 1) + d / 15);
      if x < -33 then f := -x
      else if abs (x) < 1e-6 then f := 1
      else if x > 33 then f := 0
      else f := x / (exp (x) - 1);
      beta {cm3/s} := 1e-3 * y * d * f * sqrt (sqr (pi2Dp) + sqr (pi2Di));
   end;

{Particle-particle coagulation coefficient, Sahni approximation}
{Uses external function "Mobility"}
function coag (d1, {diameter of small neutral particle : nm}
{cm3/s}        d2, {diameter of large charged or neutral particle : nm}
                h, {extra distance : nm, typically 0.115}
             gcm3, {small particle density : g cm-3, typically 1.5-2}
                q, {large particle charge : e}
               aa, {small particle polarizability : angstrom^3}
               d0, {critical diameter of quantum rebound : nm, typical 2.5}
               T0, {extra temperature of quantum rebound: K, typically 300}
                T, {air temperature : K}
               mb  {air pressure : mb} : real) : real;
   {Polarizability in angstrom^3 is often estimated as equal to 
    the number of atoms in the cluster or as r^3 for large particles}
   var x, y, d, b, m1, m2, g, p, a, fe, pp, qq : real;
   begin
      x := 8.674e-23 * T; {2*pi*k*T : SI}
      d := 1e-9 * (d1 + d2 + 2 * h); {2 * delta : SI}
      m1 := 5.236e-25 * gcm3 * d1 * d1 * d1; {SI}
      m2 := 5.236e-25 * gcm3 * d2 * d2 * d2; {SI}
      b := 1e15 * {B1 + B2 : SI}
           (mobility (28.96, 0.00171, 0.3036, 44, 0.8, mb, T, gcm3, 0, d1) +
            mobility (28.96, 0.00171, 0.3036, 44, 0.8, mb, T, gcm3, 0, d2));
      g := 2 * (b / d) * sqrt (x * m1 * m2 / (m1 + m2)); {gamma = Kc / Kk}
      if d0 = 0 then p := 1 else begin
         y := (273 / (T + T0)) * exp (3 * ln (d0 / d2));
         p := sqr (y) * exp (y) / sqr (exp (y) - 1); {sticking probability}
      end;
      if (q = 0) or (aa = 0) then fe := 1 else begin
         a := aa * 1e-30; {polarizability : m^3}
         y := 119700 * T * sqr (sqr (d / 2)) / (a * sqr (q));
           {y = T* = kT/Upol, 119700 = 8 * pi * eps0 * k / e^2}
         if y > 1 then fe := 1 + 0.106 / y + 0.263 / exp ((4/3) * ln (y))
         else begin pp := sqrt (y); qq := sqrt (pp);
           fe := 1.4691 / pp - 0.341 / qq + 0.182 * y * qq + 0.059;
         end;
      end;
      coag {cm3/s} := 1e6 *  fe * x * b * d /
         (1 + g - 0.299 * g / (exp (1.1 * ln (g)) + 0.64) + g * (1 - p) / p);
   end;

{Growth rate factor GR/GRo for the first substance,
 designed as modification of function coag at non-evaporating condensation.
 GRo is the plain Knudsen growth rate explained in (Tammet, Kulmala, 2005)}
function growthfactor1 {dimensionless, uses external function "Mobility"}
              (d1, {diameter of small neutral particle : nm}
               d2, {diameter of large charged or neutral particle : nm}
                h, {extra distance : nm, typically 0.115}
             gcm3, {small particle density : g cm-3, typically 1.5-2}
                q, {large particle charge : e}
               aa, {small particle polarizability : angstrom^3}
               d0, {critical diameter of quantum rebound : nm, typical 2.5}
               T0, {extra temperature of quantum rebound: K, typically 300}
                T, {air temperature : K}
               mb, {air pressure : mb}
              yua, {first dipole enhancement coefficient}
              yub  {second dipole enhancement coefficient} : real) : real;
    {Two alternative methods can be used (don't use them simultaneously!):
     1. In case of the method of effective polarizability, the actual value
        of aa should be presented and yua = yub = 0.
        Polarizability aa in angstrom^3 is often estimated as equal to the
        number of atoms in a cluster or as equal to r^3 for large particles.
        Polarizability of a molecule of sulphuric acid is extra high:
        about 149 angstrom^3.
     2. In case of the method by Nadykto and Yu the parameter aa must be zero,
        yua = sqr(f1 - 1) / (f2 - 1) and yub = ln((f1 - 1) / (f2 - 1)),
        where f1 is the Nadykto-Yu dipole enhanchement factor for d = 1 nm
        and f2 is the Nadykto-Yu dipole enhanchement factor for d = 2 nm.
        Nadykto and Yu expected f1 = 4.35 and f2 = 1.8 for sulphuric acid
        at temperature 298 K, in this case yua = 14 and yub = 1.43}

   var x, y, d, b, m1, m2, g, p, a, fe, pp, qq, k, u : real;
   begin
      x := 8.674e-23 * T; {2*pi*k*T : SI}
      d := 1e-9 * (d1 + d2 + 2 * h); {2 * delta : SI}
      m1 := 5.236e-25 * gcm3 * d1 * d1 * d1; {SI}
      m2 := 5.236e-25 * gcm3 * d2 * d2 * d2; {SI}
      b := 1e15 * {B1 + B2 : SI}
           (mobility (28.96, 0.00171, 0.3036, 44, 0.8, mb, T, gcm3, 0, d1) +
            mobility (28.96, 0.00171, 0.3036, 44, 0.8, mb, T, gcm3, 0, d2));
      g := 2 * (b / d) * sqrt (x * m1 * m2 / (m1 + m2)); {gamma = Kc / Kk}
      if d0 = 0 then p := 1 else begin
         y := (273 / (T + T0)) * exp (3 * ln (d0 / d2));
         p := sqr (y) * exp (y) / sqr (exp (y) - 1); {sticking probability}
      end;
      if (q = 0) or (aa = 0) then fe := 1 else begin
         a := aa * 1e-30; {polarizability : m^3}
         y := 119700 * T * sqr (sqr (d / 2)) / (a * sqr (q));
           {y = T* = kT/Upol, 119700 = 8 * pi * eps0 * k / e^2}
         if y > 1 then fe := 1 + 0.106 / y + 0.263 / exp ((4/3) * ln (y))
         else begin pp := sqrt (y); qq := sqrt (pp);
           fe := 1.4691 / pp - 0.341 / qq + 0.182 * y * qq + 0.059;
         end;
      end;
      k {m3/s} := fe * x * b * d /
         (1 + g - 0.299 * g / (exp (1.1 * ln (g)) + 0.64) + g * (1 - p) / p);
      u := sqrt (3.516e-23 * T / m1);
      y := 1.273e18 * k / (d2 * d2 * u);
      if (yua = 0) or (q = 0) then growthfactor1 := y
      else growthfactor1 := (1 + yua * exp (-yub * d2)) * y;
   end;

{Growth rate factor GR/GRo for the second substance,
 designed as modification of function coag according to nano-Koehler model.
 GRo is the plain Knudsen growth rate explained in (Tammet, Kulmala, 2005)}
function growthfactor2 {dimensionless, uses external function "Mobility"}
              (d1, {diameter of small neutral particle : nm}
               d2, {diameter of large charged or neutral particle : nm}
                h, {extra distance : nm, typically 0.115}
             gcm3, {small particle density : g cm-3}
                q, {large particle charge : e}
               aa, {small particle polarizability : angstrom^3}
               d0, {critical diameter of nano-Koehler model : nm, about 3}
                p, {power of the nano-Koehler model, about 2}
                T, {air temperature : K}
               mb  {air pressure : mb} : real) : real;
  {Polarizability aa in angstrom^3 is often estimated as equal to the
   number of atoms in a cluster or as equal to r^3 for larger particles.
   NB: the second substance is usually an organic compound and has lower
   polarizability when compared with the sulphuric acid}
   var x, y, d, b, m1, m2, g, a, fe, pp, qq, k, kk, u : real;
   begin
      x := 8.674e-23 * T; {2*pi*k*T : SI}
      d := 1e-9 * (d1 + d2 + 2 * h); {2 * delta : SI}
      m1 := 5.236e-25 * gcm3 * d1 * d1 * d1; {SI}
      m2 := 5.236e-25 * gcm3 * d2 * d2 * d2; {SI}
      b := 1e15 * {B1 + B2 : SI}
           (mobility (28.96, 0.00171, 0.3036, 44, 0.8, mb, T, gcm3, 0, d1) +
            mobility (28.96, 0.00171, 0.3036, 44, 0.8, mb, T, gcm3, 0, d2));
      g := 2 * (b / d) * sqrt (x * m1 * m2 / (m1 + m2)); {gamma = Kc / Kk}
      if (q = 0) or (aa = 0) then fe := 1 else begin
         a := aa * 1e-30; {polarizability : m^3}
         y := 119700 * T * sqr (sqr (d / 2)) / (a * sqr (q));
           {y = T* = kT/Upol, 119700 = 8 * pi * eps0 * k / e^2}
         if y > 1 then fe := 1 + 0.106 / y + 0.263 / exp ((4/3) * ln (y))
         else begin pp := sqrt (y); qq := sqrt (pp);
           fe := 1.4691 / pp - 0.341 / qq + 0.182 * y * qq + 0.059;
         end;
      end;
      k {m3/s} := fe * x * b * d /
         (1 + g - 0.299 * g / (exp (1.1 * ln (g)) + 0.64));
      u := sqrt (3.516e-23 * T / m1);
      if d2 < d0  then kk := 0 {nano-Koehler coefficient}
      else kk := 1 - exp ( p * ln (d0 / d2));
      growthfactor2 := 1.273e18 * k * kk / (d2 * d2 * u);
   end;

function ion_needlesink ({s-1}
           Z {ion mobility : cm2 V-1 s-1},
     dneedle {mm, about 0.9 for Pinus Sylvestris},
           L {m-2, total length of needles in 1 m^3 of the canopy},
        wind {m s-1, inside of the canopy},
           T {air temperature : K},
          mb {air pressure : mb} : real) : real;
   var nyy, Re, Sc, kTB : real;
   begin
      kTB := 8.62e-9 * T * Z;
      nyy := 5.5e-7 * exp (1.8 * ln (T)) / mb;
      Re := (wind * dneedle / nyy) / 2000;
      Sc := nyy / kTB;
      ion_needlesink := kTB * L * (0.3 +
         0.62 * sqrt (Re) * exp (0.33333 * ln (Sc)) *
         exp (0.8 * ln (1 + exp (0.625 * ln (Re / 282000)) )) /
         sqrt (sqrt (1 + exp (0.66667 * ln (0.4 / Sc)) ))       );
   end;

function particle_needlesink ({s-1}
   dparticle {nm},
        gcm3 {particle density : g cm-3},
     dneedle {mm, about 0.9 for Pinus Sylvestris},
           L {m-2, total length of needles in 1 m^3 of the canopy},
        wind {m s-1, inside of the canopy},
           T {air temperature : K},
          mb {air pressure : mb} : real) : real;
   var nyy, Re, Sc, kTB : real;
   begin
      kTB := 1.38e-8 * T * mobility (28.96, 0.00171, 0.3036, 44, 0.8,
                                     mb, T, gcm3, 0, dparticle);
      nyy := 5.5e-7 * exp (1.8 * ln (T)) / mb;
      Re := (wind * dneedle / nyy) / 2000;
      Sc := nyy / kTB;
      particle_needlesink := kTB * L * (0.3 +
         0.62 * sqrt (Re) * exp (0.33333 * ln (Sc)) *
         exp (0.8 * ln (1 + exp (0.625 * ln (Re / 282000)) )) /
         sqrt (sqrt (1 + exp (0.66667 * ln (0.4 / Sc)) ))       );
   end;

function parabol (x, y0, y1, y2, bottom : real) : real;
   {y0 = y (x=0), y1 = y (x=1/2), y2 = y (1)}
   var y : real;
   begin
   y := (2 * (y2 - 2 * y1 + y0) * x + 4 * y1 - y2 - 3 * y0) * x + y0;
   if y < bottom then parabol := bottom else parabol := y;
   end;

{============================================================================}
{Units:
   time : s
   particle diameter : nm
   concentration : cm-3 for charged and 100 cm-3 for neutral particles
   ionization and nucleation rate : cm-3 s-1
   recombination, coagulation and attachment coefficients : cm3 s-1
   electrical mobility : cm-2 V-1 s-1
   diffusion factor : nm-1 cm3 s-1
   growth rate : nm s-1
   charge : e
   diameter length : nm cm-3
   sink : s-1}
{Abbreviations: _p = positive,
                _n = negative,
                _z = zero charge,
                f_ = forest,
                d_ = diameter,
                n_ = concentration}

{GLOBAL} const intervalspace = 1440;  {technical limit for numberofintervals}
                 minutespace = 1440; {technical limit for numberofminutes}
               fractionspace = 100;  {technical limit for numberoffractions}
        minimumconcentration = 1e-4; {technical limit}
               baseofunknown = -1e6; {used only when formatting the output}
               nameofcontrol = 'burstcontrol.txt';
            numberofcontrols = 98;   {how many parameters in the control text}
                 spreadconst = 1.1;  {technical parameter of calculations}
                       legal = -1e3; {allowed negative value of the errors}
                  theendline = '=== THE END ==='; {finish of the control text}

{GLOBAL} var
         gr_file, {growth table file}
         control, {input file (control strings)}
         outfile : {output file}
                  text;
   indeterminacy, {symbol of unknown value defined in the control file}
   nameofresult : {name for the output file}
                  string;
   controlline,
   controlvalue : array [1..numberofcontrols + 1] of string;
         forest :  {considered or not}
      boolean;
          time0 : real;
{variables for local&temporary use}
      i, j,
      k, p : integer;
   a, b, c : real;
   d, g, z : real;
      x, y : real;
      r, s : string;
{counters of size and time}
   section,     {1..numberofsections, section is size interval in calculations}
   fraction,    {1..numberoffractions, fraction is size interval in the output}
   moment,      {-#..0..numberofmoments, moment is time step in calculations}
   minute,      {1..numberofminutes}
   interval,    {1..numberofintervals, interval is time step in presentation}
   process:     {0..numberofintervals, processes are used in forest processing}
      integer;
{parameters of size-time structure}
   numberofsections,
   numberofintervals,
   stepofprogress,
   numberofminutes,
   numberofmoments,
   momentsinminute,
   momentsininterval,
   minutesininterval : integer;
   acceleration, {coefficient, which increases the section width, usually 1}
   dd, dt : real;
{parameters of presentation structure}
   numberoffractions,
   diameterormobility,
   lin1log2,
   distr_type,
   delimiter : integer;
   d_min, d_max, {the presentation limits do not affect the simulation process}
   d_low, d_high,
   d_down, d_up,
   timeunit : real;
   d_fr,
   d_ctr,
   z_fr,
   z_ctr,
   buffer : array [0..fractionspace] of real;
   deliverstandardtable,
   deliverparameterlist,
   deliverdistribution,
   deliverplottable : boolean;
{constant environmental parameters}
   celsius,
   kelvin,
   millibar : real;
{parameters of nucleationfunction}
   riseshape,
   dropshape : integer;
   riseminutes,
   flatminutes,
   dropminutes : real;
   nucleationmoments : integer;
   topnucleation_p, f_topnucleation_p,
   topnucleation_n, f_topnucleation_n,
   topnucleation_z, f_topnucleation_z : real;
{growth factors for substance 1}
   growthunitdensity,
   extradistance,
   growthunitdiameter1,
   polarizability1,
   yu_enhancement_a,
   yu_enhancement_b,
   growthrate10, f_growthrate10,
   growthrate11, f_growthrate11,
   growthrate12, f_growthrate12,
   quantcriticalsize,
   extratemperature : real;
   n_growthtable : integer;
   growthtable : array [1..99, 1..2] of real;
{growth factors for substance 2}
   growthunitdiameter2,
   polarizability2,
   growthrate20, f_growthrate20,
   growthrate21, f_growthrate21,
   growthrate22, f_growthrate22,
   koehlercriticalsize,
   koehlerpower : real;
{parameters of ionization, background aerosol, and forest}
   ionizationrate0, f_ionizationrate0,
   ionizationrate1, f_ionizationrate1,
   ionizationrate2, f_ionizationrate2,
   backgrounddiameter0, backgroundconcentration0,
   backgrounddiameter1, backgroundconcentration1,
   backgrounddiameter2, backgroundconcentration2,
   wind0, wind1, wind2,
   pathofairinforest,
   needlediameter, needlelength : real;
{controls of the ion and particle evolution process}
   birthsize, finalsize,
   ionmobility_p, ionmobility_n,
   iondensity,
   recombinationcoefficient,
   d_ion_p, d_ion_n,
   eeta_p, eeta_n,
   chargedia,
   effectivedialength, correctionfactor,
   g_sigma_p, g_sigma_n : real;
   nucleationratio : array of real;
   n_ion_p, n_ion_n,
   ionizationrate, f_ionizationrate,
   wind,
   backgrounddiameter,
   backgroundconcentration,
   ionforestsink_p, ionforestsink_n,
   minsize, maxsize : array [0..intervalspace] of real;
   residence,
   linkmoment : array [0..intervalspace] of integer;
   beta_pz, beta_nz,
   beta_pn, beta_np : array of real;
   n_particle_p,
   n_particle_n,
   n_particle_z,
   particleaerosolsink,
   particleforestsink,
   growth_q, growth_z : array of array of single;
{output data}
   np1, np2,
   dp1, dp2,
   nn1, nn2,
   dn1, dn2,
   nz1, nz2,
   dz1, dz2,
   balance1, balance2,
   backcharge : array [0..intervalspace] of single;
   table_p,
   table_n,
   table_z : array [0..fractionspace, 0..intervalspace] of single;
   header_p, header_n, header_z : string;
{output format controls}
   ranknucl_p, ranknucl_n, ranknucl_z,
   f_ranknucl_p, f_ranknucl_n, f_ranknucl_z,
   rankion_p, rankion_n,
   rankd_p, rankd_p2,
   rankn_p, rankn_p2,
   rankd_n, rankd_n2,
   rankn_n, rankn_n2,
   rankd_z, rankd_z2,
   rankn_z, rankn_z2,
   ranktable_p, ranktable_n, ranktable_z,
   rankbalance1, rankbalance2,
   rank_nd, rank_nq, rank_q,
   ionunit, chargedunit, neutralunit,
   backchargeunit, backgroundunit,
   diatablecomposition : integer;
{LOCAL} var {used in evolution calculation}
   first, last,
   maxlast : integer;
   lowsize, topsize,
   dip, nip, nip0,
   din, nin, nin0,
   nucleation_p, nucleation_n, nucleation_z,
   ionaerosolsink_p, ionaerosolsink_n : real;
   trn, npp, npn, npz : array of real;
   dbackcharge, backgroundcharge,
   back0 : real;

function t_value (i : integer) : string;
   var k : integer;
   begin
      k := pos (' ', controlline [i]) - 1;
      if k > 0 then t_value := copy (controlline [i], 1, k)
               else t_value := '';
   end;

function i_value (i : integer) : integer;
   var k, c, y : integer;
   begin
      k := pos (' ', controlline [i]) - 1;
      if k > 0 then val (copy (controlline [i], 1, k), y, c)
               else begin y := 0; c := 1 end;
      i_value := y;
      if c <> 0 then begin
         writeln;
         writeln ('Control line ', i, ' is defective:');
         writeln ('"', controlline [i], '"');
         writeln ('Press ENTER for escape!');
         readln; halt;
      end;
   end;

function d_value (i : integer) : real;
   var k, c : integer;
          y : real;
   begin
      k := pos (' ', controlline [i]) - 1;
      if k > 0 then val (copy (controlline [i], 1, k), y, c)
               else begin y := 0; c := 1 end;
      d_value := y;
      if c <> 0 then begin
         writeln;
         writeln ('Control line ', i, ' is defective:');
         writeln ('"', controlline [i], '"');
         writeln ('Press ENTER for escape!');
         readln; halt;
      end;
   end;

procedure d_dualvalue (i : integer; var a, b : real);
   var v, k, c : integer;
   begin
      c := 1;
      v := pos ('/', controlline [i]);
      k := pos (' ', controlline [i]) - 1;
      if k > 0 then begin
         if (v > 1) and ((k - v) > 0) then begin
            val (copy (controlline [i], 1, v - 1), a, c);
            if c = 0 then val (copy (controlline [i], v + 1, k - v), b, c);
         end
         else begin
            val (copy (controlline [i], 1, k), a, c); b := a;
         end;
      end;
      if c <> 0 then begin
         writeln;
         writeln ('Control line ', i, ' is defective:');
         writeln ('"', controlline [i], '"');
         writeln ('Press ENTER for escape!');
         readln; halt;
      end;
   end;

procedure writeresult1 (x : real);
   begin
      if x > (baseofunknown + 1e3) then write (outfile, #09, x:3:1)
      else if indeterminacy = '-' then write (outfile, #09)
      else if indeterminacy = '#' then write (outfile, #09,
                                               (baseofunknown - x):3:1)
      else write (outfile, #09, indeterminacy);
   end;

procedure writeresult2 (x : real);
   begin
      if x > (baseofunknown + 1e3) then write (outfile, #09, x:4:2)
      else if indeterminacy = '-' then write (outfile, #09)
      else if indeterminacy = '#' then write (outfile, #09,
                                              (baseofunknown - x):4:2)
      else write (outfile, #09, indeterminacy);
   end;

procedure writeadjust (i : integer);
   var       b : integer; // index of fraction upper border
       n, p, q : real;
   begin
   if diameterormobility = 1 then b := i else b := numberoffractions + 1 - i;
   n := buffer [b]; // fraction concentration
   if diameterormobility = 1 then begin p := d_fr [b - 1]; q := d_fr [b] end
                             else begin q := z_fr [b - 1]; p := z_fr [b] end;
   if distr_type = 2 then n := n / (q - p);
   if distr_type = 3 then n := n / (0.4343 * ln (q / p));
   writeresult1 (n);
   end;

procedure failure (code : integer);
   begin
   writeln;
   writeln ('Sorry, this program is not able to process some specific tasks,');
   writeln ('especially when the ionization rate is very high and the number');
   writeln ('of evolution steps in a minute is low.');
   writeln;
   writeln ('You can help to develop the program mailing the failure info');
   writeln ('and a copy of the control file to hannes.tammet@ut.ee.');
   writeln;
   writeln ('The info is: Failure ', code, ' in ', revision, '.');
   writeln;
   writeln ('Press ENTER for escape!');
   readln;
   halt;
   end;

function diameter1 (z : real) : real;
   begin
   diameter1 := massdiameter (millibar, kelvin, iondensity, 1, 0.624 * z);
   end;

function tablegrowthrate (diameter : real) : real;
   var i, k : integer;
       a, b : real;
   begin
   if diameter <= growthtable [1, 1]
      then tablegrowthrate := growthtable [1, 2]
   else if diameter >= growthtable [n_growthtable, 1]
      then tablegrowthrate := growthtable [n_growthtable, 2]
   else begin
      k := 1;
      for i := 1 to n_growthtable - 1 do
      if growthtable [i, 1] < diameter then k := i;
      // diameter is between growthtable [k, 1] and growthtable [k + 1, 1]
      // linear interpolation
      a := diameter - growthtable [k, 1];
      b := growthtable [k + 1, 1] - diameter;
      tablegrowthrate :=
      (b * growthtable [k, 2] + a * growthtable [k + 1, 2]) / (a + b);
   end;
   end;

function freeairgrowthrate (diameter, time, charge : real) : real;
   begin
   if n_growthtable > 1
   then freeairgrowthrate := tablegrowthrate (diameter)
   else freeairgrowthrate :=
      parabol (time, growthrate10, growthrate11, growthrate12, 0) *
      growthfactor1 (growthunitdiameter1, diameter, extradistance,
                     growthunitdensity, charge, polarizability1,
                     quantcriticalsize, extratemperature, kelvin, millibar,
                     yu_enhancement_a, yu_enhancement_b) +
      parabol (time, growthrate20, growthrate21, growthrate22, 0) *
      growthfactor2 (growthunitdiameter2, diameter, extradistance,
                     growthunitdensity, charge, polarizability2,
                     koehlercriticalsize, koehlerpower, kelvin, millibar);
   end;

function forestgrowthrate (diameter, time, charge : real) : real;
   begin
   if n_growthtable > 1
   then forestgrowthrate := tablegrowthrate (diameter)
   else forestgrowthrate :=
      parabol (time, f_growthrate10, f_growthrate11, f_growthrate12, 0) *
      growthfactor1 (growthunitdiameter1, diameter, extradistance,
                     growthunitdensity, charge, polarizability1,
                     quantcriticalsize, extratemperature, kelvin, millibar,
                     yu_enhancement_a, yu_enhancement_b) +
      parabol (time, f_growthrate20, f_growthrate21, f_growthrate22, 0) *
      growthfactor2 (growthunitdiameter2, diameter, extradistance,
                     growthunitdensity, charge, polarizability2,
                     koehlercriticalsize, koehlerpower, kelvin, millibar);
   end;

Procedure check_file (c : char);
   var ext, s, x : string;
               f : text;
   begin
   if c = 'L' then ext := '.txt' else ext := '.xl';
   s := c + copy (nameofresult, 2, 3) + ext;
   if fileexists (s) then begin
      writeln;
      if fileisreadonly (s) then begin
         writeln ('File ', s, ' exists and is read-only');
         writeln ('Press ENTER for escape and try again!');
         readln; halt;
      end
      else begin
         {$I-} assignfile (f, s); append (f); {$I+}
         if IOresult <> 0 then begin
            writeln;
            writeln ('File ', s, ' is locked for editing by another program');
            writeln ('Press ENTER for escape and try again!');
            readln; halt;
         end;
         closefile (f);
         write (s, ' exists, owerwrite (y/n)? ');
         readln (x);
         if (x<> '') and (x<> 'y') and (x<> 'Y') then begin
            writeln ('Press ENTER for escape and try again!');
            readln; halt;
         end;
      end;
   end;
   end;

{============================================================================}

begin
try
   {for i := 150 to 250 do writeln (i, chr (i));}
   if lastyear = 9999 then writeln ('Burstsimulator ', revision)
   else begin
      writeln ('Burstsimulator ', revision, ', valid ',
               (encodedate (lastyear, lastmonth, lastday) - date):1:0, ' days');
      if date > encodedate (lastyear, lastmonth, lastday) then begin
         writeln;
         writeln ('   Sorry, the validity of the program is expired.');
         writeln;
         writeln ('   Press ENTER to escape!');
         readln; halt;
      end
   end;

{READ and analyze the control file}
   if paramcount = 0 then s := nameofcontrol
   else begin
      s := paramstr (1); k := length (s);
      repeat k := k - 1 until (s [k] = '\') or (k < 3);
      if k < 3 then begin
         writeln;
         writeln ('Wrong item was dragged, press ENTER for escape!');
         readln; halt;
      end;
      writeln ('Control file: ', s);
      chdir (copy (s, 1, k - 1));
   end;
   assignfile (control, s);
   {$I-} reset (control); {$I+}
   if IOresult <> 0 then begin
      writeln;
      writeln ('Control file is not available, press ENTER for escape!');
      readln; halt;
   end;
   k := 0;
   {$I-}
   repeat
      readln (control, s);
      if IOresult <> 0 then begin
         writeln;
         writeln ('Hardware error reading the control file,',
                  ' press ENTER for escape!');
         readln; halt;
      end;
      s := trimleft (s);
      if length (s) > 1 then begin
         if (s [1] <> '=') then begin
            k := k + 1;
            controlline [k] := s;
         end;
      end;
   until (k > numberofcontrols) or (s = theendline) or eof (control);
   {$I+}
   closefile  (control);
   if s <> theendline then begin
      writeln;
      writeln ('CORRUPTED CONTROL FILE!');
      writeln ('Press ENTER for escape!');
      readln; halt;
   end;
   if k <> numberofcontrols then begin
      writeln;
      writeln ('Wrong number of control lines in the control file!');
      writeln ('Press ENTER for escape!');
      readln; halt;
   end;

{READ the growthtable}
   n_growthtable := 0; // mark that the table is not presented
   if fileexists ('burstgrowthtable.txt') then begin
      assignfile (gr_file, 'burstgrowthtable.txt'); reset (gr_file);
      n_growthtable := 0;
      while not eof (gr_file) do begin
         n_growthtable := n_growthtable + 1;
         readln (gr_file, growthtable [n_growthtable, 1], x);
         growthtable [n_growthtable, 2] := x / 3600;
         if (growthtable [n_growthtable, 1] = 0) and (x = 0)
         then n_growthtable := n_growthtable - 1; // cancel empty lines
      end;
      closefile (gr_file);
      if n_growthtable < 2 then begin
         writeln;
         writeln ('burstgrowthtable too short (should be at least 4 lines)!');
         writeln ('Press ENTER for escape!');
         readln; halt;
      end;
      writeln ('burstgrowthtable loaded');
   end;

{PICK VALUES of control parameters}
   numberofminutes := i_value (1);
   momentsinminute := i_value (2);
   acceleration := d_value (3);
   celsius := d_value (4); kelvin := 273.15 + celsius;
   millibar := d_value (5);
   d_dualvalue (6, ionizationrate0, f_ionizationrate0);
   d_dualvalue (7, ionizationrate1, f_ionizationrate1);
   d_dualvalue (8, ionizationrate2, f_ionizationrate2);
   ionmobility_p := d_value (9);
   ionmobility_n := d_value (10);
   recombinationcoefficient := d_value (11);
   iondensity := d_value (12);
   birthsize := d_value (13);
   d_dualvalue (14, topnucleation_p, f_topnucleation_p);
   d_dualvalue (15, topnucleation_n, f_topnucleation_n);
   d_dualvalue (16, topnucleation_z, f_topnucleation_z);
   riseminutes := d_value (17);
   riseshape := i_value (18);
   flatminutes := d_value (19);
   dropminutes := d_value (20);
   dropshape := i_value (21);
   growthunitdensity := d_value (22);
   extradistance := d_value (23);
   growthunitdiameter1 := d_value (24);
   polarizability1 := d_value (25);
   a := d_value (26);
   b := d_value (27);
   if a <= 1 then begin yu_enhancement_a := 0; yu_enhancement_b := 0 end
   else begin
      if (a <= b) or (b < 1) then begin
         writeln ('Nadykto-Yu parameters do not fit the requirements');
         writeln ('Press ENTER!');
         readln; halt;
      end;
      yu_enhancement_a := sqr(a - 1) / (b - 1);
      yu_enhancement_b := ln((a - 1) / (b - 1));
   end;
   d_dualvalue (28, x, y);
   growthrate10 := x / 3600;
   f_growthrate10 := y / 3600;
   d_dualvalue (29, x, y); growthrate11 := x / 3600; f_growthrate11 := y / 3600;
   d_dualvalue (30, x, y); growthrate12 := x / 3600; f_growthrate12 := y / 3600;
   quantcriticalsize := d_value (31);
   extratemperature := d_value (32);
   growthunitdiameter2 := d_value (33);
   polarizability2 := d_value (34);
   d_dualvalue (35, x, y); growthrate20 := x / 3600; f_growthrate20 := y / 3600;
   d_dualvalue (36, x, y); growthrate21 := x / 3600; f_growthrate21 := y / 3600;
   d_dualvalue (37, x, y); growthrate22 := x / 3600; f_growthrate22 := y / 3600;
   koehlercriticalsize := d_value (38);
   koehlerpower := d_value (39);
   backgrounddiameter0 := d_value (40);
   backgrounddiameter1 := d_value (41);
   backgrounddiameter2 := d_value (42);
   backgroundconcentration0 := d_value (43);
   backgroundconcentration1 := d_value (44);
   backgroundconcentration2 := d_value (45);
   pathofairinforest := d_value (46);
   wind0 := d_value (47);
   wind1 := d_value (48);
   wind2 := d_value (49);
   needlediameter := d_value (50);
   needlelength := d_value (51);
   forest := (pathofairinforest > 0) and (needlelength > 0);
   if not forest then pathofairinforest := 0;
   if not forest then needlediameter := 1; // avoid x / 0
   diameterormobility := i_value (52);
   d_min := d_value (53);
   if diameterormobility = 2 then d_min := diameter1 (d_min);
   d_max := d_value (54);
   if diameterormobility = 2 then d_max := diameter1 (d_max);
   if d_min > d_max then begin x := d_min; d_min := d_max; d_max := x end;
   numberoffractions := i_value (55);
   if numberoffractions < 1 then numberoffractions := 1;
   lin1log2 := i_value (56);
   distr_type := i_value (57);
   minutesininterval := i_value (58);
   timeunit := i_value (59); {expressed in minutes}
   delimiter := i_value (60);
   s := uppercase (t_value (61));
   deliverstandardtable := pos ('T', s) > 0;
   deliverparameterlist := pos ('L', s) > 0;
   deliverdistribution := pos ('D', s) > 0;
   deliverplottable := pos ('P', s) > 0;
   if not (deliverstandardtable or deliverparameterlist
        or deliverdistribution or deliverplottable) then begin
      writeln ('Nothing to deliver. Press ENTER for escape!');
      readln; halt;
   end;
   ranknucl_p := i_value (62);
   ranknucl_n := i_value (63);
   ranknucl_z := i_value (64);
   f_ranknucl_p := i_value (65);
   f_ranknucl_n := i_value (66);
   f_ranknucl_z := i_value (67);
   rankion_p := i_value (68);
   rankion_n := i_value (69);
   rankd_p := i_value (70);
   rankn_p := i_value (71);
   rankd_n := i_value (72);
   rankn_n := i_value (73);
   rankd_z := i_value (74);
   rankn_z := i_value (75);
   rankbalance1 := i_value (76);
   rankd_p2 := i_value (77);
   rankn_p2 := i_value (78);
   rankd_n2 := i_value (79);
   rankn_n2 := i_value (80);
   rankd_z2 := i_value (81);
   rankn_z2 := i_value (82);
   rankbalance2 := i_value (83);
   rank_q := i_value (84);
   rank_nq := i_value (85);
   rank_nd := i_value (86);
   ranktable_p := i_value (87);
   ranktable_n := i_value (88);
   ranktable_z := i_value (89);
   ionunit := i_value (90);
   chargedunit := i_value (91);
   neutralunit := i_value (92);
   backchargeunit := i_value (93);
   backgroundunit := i_value (94);
   g_sigma_p := d_value (95);
   g_sigma_n := d_value (96);
   diatablecomposition := i_value (97);
   indeterminacy := t_value (98);
   if numberoffractions > fractionspace then begin
      writeln ('Too many size fractions, the limit is ', fractionspace, '.');
      writeln ('Press ENTER!');
      readln; halt;
   end;
   if acceleration < 1 then begin
      writeln ('Acceleration must be at least 1 in the control file!');
      writeln ('Press ENTER!');
      readln; halt;
   end;
// writeln (3600*freeairgrowthrate (3, 0.5, 0):0:3);

{ANALYZE number and name of output files}
   repeat
      write ('Please write number 1..999 for output files: ');
      readln (k);
   until (k > 0) and (k < 1000);
   str (1000 + k, nameofresult);
   if deliverstandardtable then check_file ('T');
   if deliverparameterlist then check_file ('L');
   if deliverdistribution then check_file ('D');
   if deliverplottable then check_file ('P');

{DETACH CONTROL VALUES AND IDENTIFIERS}
   for i := 1 to numberofcontrols do begin
      k := pos (' ', controlline [i]);
      if k > 1 then controlvalue [i] := copy (controlline [i], 1, k - 1);
   end;
   for i := 1 to numberofcontrols do begin
      s := controlline [i];
      if length (s) > 10 then begin
         s := copy (s, length (s) - 10, 11);
         k := pos (':', s);
         if k > 0 then controlline [i] := copy (s, k + 1, 10)
                  else controlline [i] := '';
      end;
   end;
   {as a result the control lines consist of identifiers only}

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{DETERMINE and check parameters of time structure}
   numberofintervals := numberofminutes div minutesininterval;
   stepofprogress := numberofintervals div 75 + 1;
   numberofminutes := numberofintervals * minutesininterval; {correction!}
   numberofmoments := numberofminutes * momentsinminute;
   momentsininterval := momentsinminute * minutesininterval;
   dt := 60 / momentsinminute;
   nucleationmoments := round ((riseminutes + flatminutes + dropminutes)
                       * momentsinminute);
   if numberofintervals > intervalspace then begin
      writeln ('Too many intervals, the limit is ', intervalspace, '.');
      writeln ('Press ENTER!');
      readln; halt;
   end;

{DETERMINE section width dd and numberofsections}
   finalsize := birthsize; {initial value of the increasing particle size}
   dd := 0;
   a := numberofmoments;
   if forest then a := a - pathofairinforest / (wind2 * dt);
   for moment := 1 to numberofmoments do begin
      b := freeairgrowthrate (finalsize, moment / numberofmoments, 1) * dt;
      c := b;
      if b > dd then dd := b;
      if forest then begin
         c :=  forestgrowthrate (finalsize, moment / numberofmoments, 1) * dt;
         if c > dd then dd := c;
      end;
      if moment > a then finalsize := finalsize + c
                    else finalsize := finalsize + b;
   end;
   dd := acceleration * dd;
   finalsize := birthsize + 1.05 * (finalsize - birthsize); {spread reserve}
   numberofsections := 1 + round ((finalsize - birthsize) / dd);
   {maximum diameter in calculations is finalsize, the size d_max
    may be different, it is used in fraction presentation only}
   if numberofsections < 10 then begin
      writeln ('Too low number of sections (', numberofsections,
               '), revise the control file!');
      writeln ('Press ENTER!');
      readln; halt;
   end;

{INFO}
   writeln;
   writeln (numberofmoments, ' time steps and ',
            numberofsections, ' size sections up to ',
            finalsize:4:2, ' nm');
   for interval := 0 to numberofintervals do
      if interval mod stepofprogress = 0 then write ('_');
   writeln;
   time0 := now;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}
{INITIALIZATION of dynamic arrays}
   setlength (n_particle_p, numberofsections + 1, numberofintervals + 1);
   setlength (n_particle_n, numberofsections + 1, numberofintervals + 1);
   setlength (n_particle_z, numberofsections + 1, numberofintervals + 1);
   setlength (growth_q, numberofsections + 1, numberofintervals + 1);
   setlength (growth_z, numberofsections + 1, numberofintervals + 1);
   setlength (particleaerosolsink, numberofsections + 1, numberofintervals + 1);
   setlength (particleforestsink, numberofsections + 1, numberofintervals + 1);
   setlength (trn, numberofsections + 1);
   setlength (npp, numberofsections + 1);
   setlength (npn, numberofsections + 1);
   setlength (npz, numberofsections + 1);
   setlength (beta_pz, numberofsections + 1);
   setlength (beta_nz, numberofsections + 1);
   setlength (beta_pn, numberofsections + 1);
   setlength (beta_np, numberofsections + 1);
   setlength (nucleationratio, numberofmoments + 1);

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}
{INITIALIZATION of constant variables}
   chargedia := 1.671E4 / kelvin; {e^2/(4*pi*eps0*k*T) : nm}
   d_ion_p :=  diameter1 (ionmobility_p); {nm}
   d_ion_n :=  diameter1 (ionmobility_n); {nm}
   eeta_p := 5.414E-11 * kelvin * ionmobility_p; {2*pi*D=(2*pi*k/e)*T*Z}
   eeta_n := 5.414E-11 * kelvin * ionmobility_n; {nm-1 cm3 s-1}

{PRECALCULATION of nucleationratio}
   for moment := 0 to numberofmoments do begin
      x := moment / momentsinminute; {x = continuous minute}
      if x < riseminutes then begin
         y := x / riseminutes;
         if riseshape > 1 then y := sin (1.5708 * y);
         if riseshape = 3 then y := sqr (y);
      end
      else if x < (riseminutes + flatminutes) then y := 1
      else if moment < nucleationmoments then begin
         if dropminutes = 0 then y := 0
         else y := (nucleationmoments - moment) /
                   (dropminutes * momentsinminute);
         if dropshape > 1 then y := sin (1.5708 * y);
         if dropshape = 3 then y := sqr (y);
      end
      else y := 0;
      nucleationratio [moment] := y;
   end;

{PRECALCULATION of wind, air residence time in forest,
 background aerosol, and ionization rates in free air and in forest}
   for interval := 0 to numberofintervals do begin
      x := (interval + 0.5) / numberofintervals;
      wind [interval] :=
      parabol (x, wind0, wind1, wind2, 0.1);
      residence [interval] :=
         round (pathofairinforest / (wind [interval] * dt));
      backgrounddiameter [interval] :=
         parabol (x,
         backgrounddiameter0, backgrounddiameter1, backgrounddiameter2, 20);
      backgroundconcentration [interval] :=
         parabol (x, backgroundconcentration0, backgroundconcentration1,
         backgroundconcentration2, 0);
      ionizationrate [interval] :=
         parabol (x, ionizationrate0, ionizationrate1, ionizationrate2, 0);
      f_ionizationrate [interval] :=
         parabol (x, f_ionizationrate0,
         f_ionizationrate1, f_ionizationrate2, 0);
   end;

{PRECALCULATION of linkmoments}
   for process := 0 to numberofintervals do linkmoment [process] :=
      process * momentsininterval - residence [process];

{PRECALCULATION of ion sinks on forest and particle sinks on background aerosol}
   for interval := 0 to numberofintervals do begin
      if interval mod stepofprogress = 0 then write (#176);
      ionforestsink_p [interval] := ion_needlesink (ionmobility_p,
            needlediameter, needlelength, wind [interval], kelvin, millibar);
      ionforestsink_n [interval] := ion_needlesink (ionmobility_n,
            needlediameter, needlelength, wind [interval], kelvin, millibar);
      for section := 1 to numberofsections do
         particleaerosolsink [section, interval] :=
            coag (birthsize + section * dd,
            backgrounddiameter [interval], extradistance,
            growthunitdensity, 0, 0, quantcriticalsize, 0,
            kelvin, millibar) * backgroundconcentration [interval];
   end;
   writeln;

{PRECALCULATION of attachment coefficients : cm3 s-1}
{beta_pz, beta_nz, beta_pn, beta_np}
   for section := 1 to numberofsections do begin
      d := birthsize + section * dd;
      beta_pz [section] := beta ( 0, d_ion_p, iondensity, d, kelvin, millibar);
      beta_nz [section] := beta ( 0, d_ion_n, iondensity, d, kelvin, millibar);
      beta_pn [section] := beta (-1, d_ion_p, iondensity, d, kelvin, millibar);
      beta_np [section] := beta (-1, d_ion_n, iondensity, d, kelvin, millibar);
   end;

{PRECALCULATION of free air growth during 1 second
 divided to the section width dd}
   write (#177); {growth rate in interval < 1 is not used}
   for interval := 1 to numberofintervals do begin
      if interval mod stepofprogress = 0 then write (#177);
      x := (interval + 0.5) / numberofintervals;
      for section := 1 to numberofsections do begin
         d := birthsize + section * dd;
         growth_q [section, interval] := freeairgrowthrate (d, x, 1) / dd;
         growth_z [section, interval] := freeairgrowthrate (d, x, 0) / dd;
      end;
   end;
   writeln;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{FREE AIR EQUILIBRIUM ion concentrations and background aerosol charge}
  {first approximation}
   backgroundcharge := 0;
   effectivedialength :=
        (backgrounddiameter [0] - 1.5) * backgroundconcentration [0];
   correctionfactor := (backgrounddiameter [0] + 9) /
                       (backgrounddiameter [0] + 23);
   a := recombinationcoefficient;
   b := 0.5 * (eeta_p + eeta_n) * effectivedialength;
   c := ionizationrate [0];
   if a = 0 then nip := c / b
   else nip := (sqrt (b * b + 4 * a * c) - b) / (2 * a);
   nin := nip;
  {iteration}
   c := 10; {time step during iteration}
   repeat
      ionaerosolsink_p := eeta_p * (effectivedialength -
         correctionfactor * chargedia * backgroundcharge);
      ionaerosolsink_n := eeta_n * (effectivedialength +
         correctionfactor * chargedia * backgroundcharge);
      dbackcharge := ionaerosolsink_p * nip - ionaerosolsink_n * nin;
      x := ionizationrate [0] - recombinationcoefficient * nip * nin;
      dip := x - ionaerosolsink_p * nip;
      din := x - ionaerosolsink_n * nin;
      backgroundcharge := backgroundcharge + c * dbackcharge;
      nip := nip + c * dip;
      nin := nin + c * din;
   until (abs (dip) < 1e-6)
     and (abs (din) < 1e-6)
     and (abs (dbackcharge) < 1e-6);

{PRENUCLEATION STATUS in FREE AIR}
   k := 0;
   while linkmoment [k + 1] <= 0 do k := k + 1;
   for process := 0 to k do begin;
      write (process mod 10);
      n_ion_p [process] := nip;
      n_ion_n [process] := nin;
      backcharge [process] := backgroundcharge;
      minsize [process] := 1;
      maxsize [process] := 10; {10 sections of spread reserve}
      for section := 1 to numberofsections do begin
         n_particle_p [section, process] := 0;
         n_particle_n [section, process] := 0;
         n_particle_z [section, process] := 0;
      end;
   end;
   for section := 1 to numberofsections do begin
      npp [section] := 0;
      npn [section] := 0;
      npz [section] := 0;
   end;

{FREE AIR EVOLUTION of IONS and PARTICLES}
   {initial values of limits of nonzero region expressed in dd,
    duration of a moment = dt}

   lowsize := 1; topsize := 10; {10 sections of spread reserve}
   for moment := 1  to linkmoment [numberofintervals] do begin

      interval := 1 + (moment - 1) div momentsininterval;
      effectivedialength := (backgrounddiameter [interval] - 1.5) *
                             backgroundconcentration [interval];
      correctionfactor := (backgrounddiameter [interval] + 9) /
                          (backgrounddiameter [interval] + 23);

     {limits of nonzero region}
      k := trunc (lowsize);
      if k > numberofsections then k := numberofsections;
      if moment > nucleationmoments then lowsize :=
                 lowsize + (growth_z [k, interval] / spreadconst) * dt;
      if lowsize > numberofsections - 1 then lowsize := numberofsections - 1;
      first := trunc (lowsize);
      k := trunc (topsize);
      if k > numberofsections then k := numberofsections;
      topsize := topsize + (growth_q [k, interval] * spreadconst) * dt;
      if topsize > numberofsections then topsize := numberofsections;
      last := trunc (topsize);
      maxlast := last;

     {rates of ion-induced nucleation that may be corrected when too large}
      nucleation_p := topnucleation_p * nucleationratio [moment];
      nucleation_n := topnucleation_n * nucleationratio [moment];
      nucleation_z := topnucleation_z * nucleationratio [moment];

     {sinks of ions on background aerosol}
      ionaerosolsink_p := eeta_p * (effectivedialength -
         correctionfactor * chargedia * backgroundcharge);
      ionaerosolsink_n := eeta_n * (effectivedialength +
         correctionfactor * chargedia * backgroundcharge);

     {evolution of charge of background aerosol}
      dbackcharge := ionaerosolsink_p * nip - ionaerosolsink_n * nin;
      for section := first to last do dbackcharge := dbackcharge +
         particleaerosolsink [section, interval] *
         (npp [section] - npn [section]);
      backgroundcharge := backgroundcharge + dbackcharge * dt;

     {di# = (increment of ion concentration) / dt}
     {evolution of positive cluster ions}
      x := ionaerosolsink_p;
      for section := first to last do
         x := x + beta_pz [section] * npz [section]
                + beta_pn [section] * npn [section];
      dip := ionizationrate [interval] - nucleation_p
           - recombinationcoefficient * nip * nin
           - x * nip;
      if (nip + dip * dt) < 0 then begin {nucleation constraint}
         dip := -nip / dt;
         nucleation_p := ionizationrate [interval] - dip
           - recombinationcoefficient * nip * nin
           - x * nip;
         if nucleation_p < legal then failure (1);
      end;
      nip := nip + dip * dt;

     {evolution of negative cluster ions}
      x := ionaerosolsink_n;
      for section := first to last do
         x := x + beta_nz [section] * npz [section]
                + beta_np [section] * npp [section];
      din := ionizationrate [interval] - nucleation_n
           - recombinationcoefficient * nip * nin
           - x * nin;
      if (nin + din * dt) < 0 then begin {nucleation constraint}
         din := -nin / dt;
         nucleation_n := ionizationrate [interval] - din
           - recombinationcoefficient * nip * nin
           - x * nin;
         if nucleation_n < legal then failure (2);
      end;
      nin := nin + din * dt;

     {evolution of positive nanometer particles}
      {1 s transfer from section}
      for section := first to last do begin
         trn [section] := growth_q [section, interval] * npp [section];
         if (trn [section] * dt) > npp [section]
         then trn [section] := (npp [section] / dt);
      end;
      if first > 1 then trn [first - 1] := trn [last]; {circular conservation}
      for section := first to last do begin
         if section = 1 then x := nucleation_p else x := trn [section - 1];
         npp [section] := npp [section] + dt * (x - trn [section]
            + beta_pz [section] * nip * npz [section]
            - (beta_np [section] * nin
               + particleaerosolsink [section, interval]) * npp [section]);
         if npp [section] < legal then failure (3);
      end; {of positive particles}

     {evolution of negative nanometer particles}
      {1 s transfer from section}
      for section := first to last do begin
         trn [section] := growth_q [section, interval] * npn [section];
         if (trn [section] * dt) > npn [section]
         then trn [section] := (npn [section] / dt);
      end;
      if first > 1 then trn [first - 1] := trn [last]; {circular conservation}
      for section := first to last do begin
         if section = 1 then x := nucleation_n else x := trn [section - 1];
         npn [section] := npn [section] + dt * (x - trn [section]
            + beta_nz [section] * nin * npz [section]
            - (beta_pn [section] * nip
               + particleaerosolsink [section, interval]) * npn [section]);
         if npn [section] < legal then failure (4);
      end; {of negative particles}

     {evolution of neutral nanometer particles}
      {1 s transfer from section}
      for section := first to last do begin
         trn [section] := growth_z [section, interval] * npz [section];
         if (trn [section] * dt) > npz [section]
         then trn [section] := (npz [section] / dt);
      end;
      if first > 1 then trn [first - 1] := trn [last]; {circular conservation}
      for section := first to last do begin
         if section = 1 then x := nucleation_z else x := trn [section - 1];
         npz [section] := npz [section] + dt * (x - trn [section]
            + beta_pn [section] * nip * npn [section]
            + beta_np [section] * nin * npp [section]
            - (beta_pz [section] * nip + beta_nz [section] * nin
               + particleaerosolsink [section, interval]) * npz [section]);
         if npz [section] < legal then failure (5);
      end; {of neutral particles}

      if first > 1 then for section := 1 to first - 1 do begin
         npp [section] := 0;
         npn [section] := 0;
         npz [section] := 0;
      end;

     {store the results at linkmoments of processes}
      for process := 1 to numberofintervals do
      if moment = linkmoment [process] then begin
         if process mod stepofprogress = 0 then write (process mod 10);
         n_ion_p [process] := nip;
         n_ion_n [process] := nin;
         backcharge [process] := backgroundcharge;
         for section := 1 to numberofsections do begin
            n_particle_p [section, process] := npp [section];
            n_particle_n [section, process] := npn [section];
            n_particle_z [section, process] := npz [section];
         end;
         minsize [process] := lowsize;
         maxsize [process] := topsize;
      end; {of store}

   end; {of free air moment counter}
   writeln;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{RECALCULATION of growth_q, growth_z [section, interval]
 and calculation of particle sinks on forest}
   if forest then write (#205); {growth rate in interval < 1 is not used}
   if forest then for interval := 1 to numberofintervals do begin
      if interval mod stepofprogress = 0 then write (#205);
      x := (interval + 0.5) / numberofintervals;
      for section := 1 to numberofsections do begin
         d := birthsize + section * dd;
         growth_q [section, interval] := forestgrowthrate (d, x, 1) / dd;
         growth_z [section, interval] := forestgrowthrate (d, x, 0) / dd;
         particleforestsink [section, interval] :=
            particle_needlesink (birthsize + section * dd, growthunitdensity,
            needlediameter, needlelength, wind [interval], kelvin, millibar);
      end;
   end;
   writeln;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{FOREST EVOLUTION of IONS and PARTICLES}
  {initial values of limits of nonzero region expressed in dd}
   if forest then
   for process := 0 to numberofintervals do begin

      if process mod stepofprogress = 0 then write (process mod 10);
     {pick interim values from the storage}
      nip := n_ion_p [process];
      nin := n_ion_n [process];
      backgroundcharge := backcharge [process];
      lowsize := minsize [process];
      topsize := maxsize [process];
      for section := 1 to numberofsections do begin
         npp [section] := n_particle_p [section, process];
         npn [section] := n_particle_n [section, process];
         npz [section] := n_particle_z [section, process];
      end;

      for moment := linkmoment [process] + 1 to process * momentsininterval do
      if moment < 1 then begin
         effectivedialength := (backgrounddiameter [0] - 1.5) *
                                backgroundconcentration [0];
         correctionfactor := (backgrounddiameter [0] + 9) /
                             (backgrounddiameter [0] + 23);
         ionaerosolsink_p := eeta_p * (effectivedialength -
            correctionfactor * chargedia * backgroundcharge);
         ionaerosolsink_n := eeta_n * (effectivedialength +
            correctionfactor * chargedia * backgroundcharge);
         dbackcharge := ionaerosolsink_p * nip - ionaerosolsink_n * nin;
         x := f_ionizationrate [0] - recombinationcoefficient * nip * nin;
         dip := x - (ionaerosolsink_p + ionforestsink_p [0]) * nip;
         din := x - (ionaerosolsink_n + ionforestsink_n [0]) * nin;
         backgroundcharge := backgroundcharge + dbackcharge * dt;
         nip := nip + dip * dt; if nip < legal then failure (6);
         nin := nin + din * dt; if nin < legal then failure (7);
         {particle concentrations remain still 0 while moment < 1}
      end
      else begin
         interval := 1 + (moment - 1) div momentsininterval;
         effectivedialength := (backgrounddiameter [interval] - 1.5) *
                                backgroundconcentration [interval];
         correctionfactor := (backgrounddiameter [interval] + 9) /
                             (backgrounddiameter [interval] + 23);

        {limits of nonzero region}
         k := trunc (lowsize);
         if k > numberofsections then k := numberofsections;
         if moment > nucleationmoments then lowsize :=
                    lowsize + (growth_z [k, interval] / spreadconst) * dt;
         if lowsize > numberofsections - 1 then lowsize := numberofsections - 1;
         first := trunc (lowsize);
         k := trunc (topsize);
         if k > numberofsections then k := numberofsections;
         topsize := topsize + (growth_q [k, interval] * spreadconst) * dt;
         if topsize > numberofsections then topsize := numberofsections;
         last := trunc (topsize);
         if last > maxlast then maxlast := last;

        {rates of ion-induced nucleation that may be corrected when too large}
         nucleation_p := f_topnucleation_p * nucleationratio [moment];
         nucleation_n := f_topnucleation_n * nucleationratio [moment];
         nucleation_z := f_topnucleation_z * nucleationratio [moment];

        {sinks of ions on background aerosol}
         ionaerosolsink_p := eeta_p * (effectivedialength -
            correctionfactor * chargedia * backgroundcharge);
         ionaerosolsink_n := eeta_n * (effectivedialength +
            correctionfactor * chargedia * backgroundcharge);

        {evolution of charge of background aerosol}
         dbackcharge := ionaerosolsink_p * nip - ionaerosolsink_n * nin;
         for section := first to last do dbackcharge := dbackcharge +
            particleaerosolsink [section, interval] *
            (npp [section] - npn [section]);
         backgroundcharge := backgroundcharge + dbackcharge * dt;

        {di# = (increment of ion concentration) / dt}
        {evolution of positive cluster ions}
         x := ionaerosolsink_p + ionforestsink_p [interval];
         for section := first to last do
            x := x + beta_pz [section] * npz [section]
                   + beta_pn [section] * npn [section];
         dip := f_ionizationrate [interval] - nucleation_p
              - recombinationcoefficient * nip * nin
              - x * nip;
         if (nip + dip * dt) < 0 then begin {nucleation constraint}
            dip := -nip / dt;
            nucleation_p := f_ionizationrate [interval] - dip
              - recombinationcoefficient * nip * nin
              - x * nip;
            if nucleation_p < legal then failure (8);
         end;
         nip := nip + dip * dt;

        {evolution of negative cluster ions}
         x := ionaerosolsink_n + ionforestsink_n [interval];
         for section := first to last do
            x := x + beta_nz [section] * npz [section]
                   + beta_np [section] * npp [section];
         din := f_ionizationrate [interval] - nucleation_n
              - recombinationcoefficient * nip * nin
              - x * nin;
         if (nin + din * dt) < 0 then begin {nucleation constraint}
            din := -nin / dt;
            nucleation_n := f_ionizationrate [interval] - din
              - recombinationcoefficient * nip * nin
              - x * nin;
            if nucleation_n < legal then failure (9);
         end;
         nin := nin + din * dt;

        {evolution of positive nanometer particles}
         {1 s transfer from section}
         for section := first to last do begin
            trn [section] := growth_q [section, interval] * npp [section];
            if (trn [section] * dt) > npp [section]
            then trn [section] := (npp [section] / dt);
         end;
         if first > 1 then trn [first - 1] := trn [last]; {circular conserv}
         for section := first to last do begin
            if section = 1 then x := nucleation_p else x := trn [section - 1];
            npp [section] := npp [section] + dt * (x - trn [section]
               + beta_pz [section] * nip * npz [section]
               - (beta_np [section] * nin
                  + particleaerosolsink [section, interval]
                  + particleforestsink [section, interval]) * npp [section]);
            if npp [section] < legal then failure (10);
         end; {of positive particles}

        {evolution of negative nanometer particles}
         {1 s transfer from section}
         for section := first to last do begin
            trn [section] := growth_q [section, interval] * npn [section];
            if (trn [section] * dt) > npn [section]
            then trn [section] := (npn [section] / dt);
         end;
         if first > 1 then trn [first - 1] := trn [last]; {circular conserv}
         for section := first to last do begin
            if section = 1 then x := nucleation_n else x := trn [section - 1];
            npn [section] := npn [section] + dt * (x - trn [section]
               + beta_nz [section] * nin * npz [section]
               - (beta_pn [section] * nip
                  + particleaerosolsink [section, interval]
                  + particleforestsink [section, interval]) * npn [section]);
            if npn [section] < legal then failure (11);
         end; {of negative particles}

        {evolution of neutral nanometer particles}
         {1 s transfer from section}
         for section := first to last do begin
            trn [section] := growth_z [section, interval] * npz [section];
            if (trn [section] * dt) > npz [section]
            then trn [section] := (npz [section] / dt);
         end;
         if first > 1 then trn [first - 1] := trn [last]; {circular conserv}
         for section := first to last do begin
            if section = 1 then x := nucleation_z else x := trn [section - 1];
            npz [section] := npz [section] + dt * (x - trn [section]
               + beta_pn [section] * nip * npn [section]
               + beta_np [section] * nin * npp [section]
               - (beta_pz [section] * nip + beta_nz [section] * nin
                  + particleaerosolsink [section, interval]
                  + particleforestsink [section, interval]) * npz [section]);
            if npz [section] < legal then failure (12);
         end; {of neutral particles}

         if first > 1 then for section := 1 to first - 1 do begin
            npp [section] := 0;
            npn [section] := 0;
            npz [section] := 0;
         end;

      end; {of moment counter}

     {save the results at the end of a process}
      n_ion_p [process] := nip;
      n_ion_n [process] := nin;
      backcharge [process] := backgroundcharge;
      for section := 1 to numberofsections do begin
         n_particle_p [section, process] := npp [section];
         n_particle_n [section, process] := npn [section];
         n_particle_z [section, process] := npz [section];
      end;

   end; {of process counter}

   if forest then writeln;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{INTERIM INFORMATION:
   birthsize, finalsize, dd, d_min, d_max
   numberofsections, numberoffractions
   diameterormobility, lin1log2, distr_type
   minutesininterval, numberofintervals, timeunit
   n_ion_p [interval], interval = 0..numberofintervals
   n_ion_n [interval], interval = 0..numberofintervals
   backcharge [interval], interval = 0..numberofintervals
   n_particle_p [section, interval], section = 1..numberofsections
   n_particle_n [section, interval], interval = 0..numberofintervals
   n_particle_z [section, interval]
   nucleationratio [minute]
UNITS:
   time : specified unit
   size : nm
   concentration : cm-3 }

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{COMPILE OUTPUT DATA}
  {all data correspond to the diameter distribution
  {borders and real size limits of limited size distribution}
   first := round ((d_min - birthsize) / dd);
   last := round ((d_max - birthsize) / dd);
   if first < 1 then first := 1;
   if last > numberofsections then last := numberofsections;
   if d_min < birthsize then d_low := birthsize else d_low := d_min;
   if d_max > finalsize then d_high := finalsize else d_high := d_max;
   d_low := baseofunknown - d_low;
   d_high := baseofunknown - d_high;
   d_down := baseofunknown - birthsize;
   d_up := baseofunknown - finalsize;
  {integral concentrations and average diameters}
   b := dd / 2;
   c := numberofintervals div 2;
   for interval := 0 to numberofintervals do begin
      x := 0; y := 0; {positive particles limited}
      for section := first to last do begin
         x := x + n_particle_p [section, interval];
         y := y + (birthsize + section * dd - b)
                  * n_particle_p [section, interval];
      end;
      if x > minimumconcentration then y := y / x
      else if interval < c then y := d_low else y := d_high;
      np2 [interval] := x;
      dp2 [interval] := y;

      x := 0; y := 0; {positive particles unlimited}
      for section := 1 to numberofsections do begin
         x := x + n_particle_p [section, interval];
         y := y + (birthsize + section * dd - b)
                  * n_particle_p [section, interval];
      end;
      if x > minimumconcentration then y := y / x
      else if interval < c then y := d_down else y := d_up;
      np1 [interval] := x;
      dp1 [interval] := y;

      x := 0; y := 0; {negative particles limited}
      for section := first to last do begin
         x := x + n_particle_n [section, interval];
         y := y + (birthsize + section * dd - b)
                  * n_particle_n [section, interval];
      end;
      if x > minimumconcentration then y := y / x
      else if interval < c then y := d_low else y := d_high;
      nn2 [interval] := x;
      dn2 [interval] := y;

      x := 0; y := 0; {negative particles unlimited}
      for section := 1 to numberofsections do begin
         x := x + n_particle_n [section, interval];
         y := y + (birthsize + section * dd - b)
                  * n_particle_n [section, interval];
      end;
      if x > minimumconcentration then y := y / x
      else if interval < c then y := d_down else y := d_up;
      nn1 [interval] := x;
      dn1 [interval] := y;
      x := 0; y := 0; {neutral particles limited}
      for section := first to last do begin
         x := x + n_particle_z [section, interval];
         y := y + (birthsize + section * dd - b)
                  * n_particle_z [section, interval];
      end;
      if x > minimumconcentration then y := y / x
      else if interval < c then y := d_low else y := d_high;
      nz2 [interval] := x;
      dz2 [interval] := y;

      x := 0; y := 0; {neutral particles unlimited}
      for section := 1 to numberofsections do begin
         x := x + n_particle_z [section, interval];
         y := y + (birthsize + section * dd - b)
                  * n_particle_z [section, interval];
      end;
      if x > minimumconcentration then y := y / x
      else if interval < c then y := d_down else y := d_up;
      nz1 [interval] := x;
      dz1 [interval] := y;
   end; {of integral concentrations and average diameters}

  {balance index of positive and negative particles}
   if (topnucleation_p + topnucleation_n) = 0
   then z := 10 * ionmobility_p / (ionmobility_p + ionmobility_n)
   else z := 10 * topnucleation_p / (topnucleation_p + topnucleation_n);
   balance1 [0] := baseofunknown - z;
   balance2 [0] := baseofunknown - z;
   for interval := 1 to numberofintervals do begin
      x := (np1 [interval] + nn1 [interval]);
      if x < 0.01 then balance1 [interval] := baseofunknown - 5
                  else balance1 [interval] := 10 * np1 [interval] / x;
      x := (np2 [interval] + nn2 [interval]);
      if x < 0.01 then balance2 [interval] := baseofunknown - 5
                  else balance2 [interval] := 10 * np2 [interval] / x;
   end; {of balance index}

  {make fraction diameter upper borders d_fr and mobility borders z_fr arranged
   according to ascendant size and considering diameterormobility and lin1log2}
   if diameterormobility = 1 then begin a := d_min; b := d_max end else begin
     {distribute mobility}
      a := 1.602 * mobility (28.96, 0.00171, 0.3036, 44, 0.8,
                             millibar, kelvin, iondensity, 1, d_min);
      b := 1.602 * mobility (28.96, 0.00171, 0.3036, 44, 0.8,
                             millibar, kelvin, iondensity, 1, d_max);
   end;
   {buffer = border for diameter of mobility}
   for fraction := 0 to numberoffractions do if lin1log2 = 1
   then buffer [fraction] := a + fraction * (b - a) / numberoffractions
   else buffer [fraction] := exp (ln (a) +
                           fraction * ln (b / a) / numberoffractions);
   for fraction := 0 to numberoffractions do
   if diameterormobility = 1 then d_fr [fraction] := buffer [fraction]
   else begin d_fr [fraction] := diameter1 (buffer [fraction]);
              z_fr [fraction] := buffer [fraction] end;
  {NB: the sections of n_particle and the borders of d as well as z
   are arranged according to the ascendant diameters}

  {make fraction centers}
   for fraction := 1 to numberoffractions do
   if lin1log2 = 1 then begin
      d_ctr [fraction] := (d_fr [fraction] + d_fr [fraction - 1]) / 2;
      z_ctr [fraction] := (z_fr [fraction] + z_fr [fraction - 1]) / 2;
   end
   else begin
      d_ctr [fraction] := sqrt (d_fr [fraction] * d_fr [fraction - 1]);
      z_ctr [fraction] := sqrt (z_fr [fraction] * z_fr [fraction - 1]);
   end;

  {make distribution subtables for fraction concentrations

    pseudoindex:       x                     y
                --+-------+-------+-------+-------+---
    section index:        i      i+1     j-1      j    }

   for interval := 0 to numberofintervals do begin
      for fraction := 1 to numberoffractions do begin {positive particles}
         x := (d_fr [fraction - 1] - birthsize) / dd; {pseudoindex}
         y := (d_fr [fraction] - birthsize) / dd; {pseudoindex}
         i := trunc (x) + 1; // from
         j := trunc (y) + 1; // until
         if (i > last) or (j < first) then c := 0 else begin
            if i < first then begin i := first; x := first end;
            if j > last then begin j := last; y := last end;
            c := (i - x) * n_particle_p [i, interval]
               + (y - j + 1) * n_particle_p [j, interval];
            if (j - 1) > (i + 1) then
            for k := i + 1 to j - 1 do
               c := c + n_particle_p [k, interval];
         end;
         table_p [fraction, interval] := c;
      end {of positive particles};

      for fraction := 1 to numberoffractions do begin {negative particles}
         x := (d_fr [fraction - 1] - birthsize) / dd; {pseudoindex}
         y := (d_fr [fraction] - birthsize) / dd; {pseudoindex}
         i := trunc (x) + 1;
         j := trunc (y) + 1;
         if (i > last) or (j < first) then c := 0 else begin
            if i < first then begin i := first; x := first end;
            if j > last then begin j := last; y := last end;
            c := (i - x) * n_particle_n [i, interval]
               + (y - j + 1) * n_particle_n [j, interval];
            if (j - 1) > (i + 1) then
            for k := i + 1 to j - 1 do
               c := c + n_particle_n [k, interval];
         end;
         table_n [fraction, interval] := c;
      end {of negative particles};

      for fraction := 1 to numberoffractions do begin {neutral particles}
         x := (d_fr [fraction - 1] - birthsize) / dd; {pseudoindex}
         y := (d_fr [fraction] - birthsize) / dd; {pseudoindex}
         i := trunc (x) + 1;
         j := trunc (y) + 1;
         if (i > last) or (j < first) then c := 0 else begin
            if i < first then begin i := first; x := first end;
            if j > last then begin j := last; y := last end;
            c := (i - x) * n_particle_z [i, interval]
               + (y - j + 1) * n_particle_z [j, interval];
            if (j - 1) > (i + 1) then
            for k := i + 1 to j - 1 do
               c := c + n_particle_z [k, interval];
         end;
         table_z [fraction, interval] := c;
      end {of neutral particles};
   end {of distribution subtables};
  {tables are still arranged according to increasing fraction size
   and include fraction concentrations}

  {make headers of distribution subtables}
   header_p := ''; header_n := ''; header_z := '';
   for fraction := 1 to numberoffractions do begin
      if distr_type = 1 then begin
         if diameterormobility = 1 then begin
            str (d_fr [fraction - 1]:4:2, s);
            str (-d_fr [fraction]:4:2, r)
         end
         else begin
            str (z_fr [numberoffractions - fraction + 1]:4:2, s);
            str (-z_fr [numberoffractions - fraction]:4:2, r)
         end;
         s := s + r;
      end
      else if diameterormobility = 1 then str (d_ctr [fraction]:4:2, s)
         else str (z_ctr [numberoffractions - fraction + 1]:4:2, s);
      if fraction < numberoffractions then s := s + #09;
      header_p := header_p + 'p' + s;
      header_n := header_n + 'n' + s;
      header_z := header_z + 'z' + s;
   end;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{WRITE OUTPUT DATA}
if deliverstandardtable then begin
   s := nameofresult; s [1] := 'T'; s := s + '.xl';
   assignfile (outfile, s); rewrite (outfile);
  {heading}
   write (outfile, 'time');
   for k := 1 to 99 do begin
      if ranknucl_p = k then write (outfile, #09, controlline [62]);
      if ranknucl_n = k then write (outfile, #09, controlline [63]);
      if ranknucl_z = k then write (outfile, #09, controlline [64]);
      if f_ranknucl_p = k then write (outfile, #09, controlline [65]);
      if f_ranknucl_n = k then write (outfile, #09, controlline [66]);
      if f_ranknucl_z = k then write (outfile, #09, controlline [67]);
      if rankion_p = k then write (outfile, #09, controlline [68],
                                           '/', ionunit);
      if rankion_n = k then write (outfile, #09, controlline [69],
                                           '/', ionunit);
      if rankd_p = k then write (outfile, #09, controlline [70]);
      if rankn_p = k then write (outfile, #09, controlline [71],
                                         '/', chargedunit);
      if rankd_n = k then write (outfile, #09, controlline [72]);
      if rankn_n = k then write (outfile, #09, controlline [73],
                                         '/', chargedunit);
      if rankd_z = k then write (outfile, #09, controlline [74]);
      if rankn_z = k then write (outfile, #09, controlline [75],
                                         '/', neutralunit);
      if rankbalance1 = k then write (outfile, #09, controlline [76]);
      if rankd_p2 = k then write (outfile, #09, controlline [77]);
      if rankn_p2 = k then write (outfile, #09, controlline [78],
                                          '/', chargedunit);
      if rankd_n2 = k then write (outfile, #09, controlline [79]);
      if rankn_n2 = k then write (outfile, #09, controlline [80],
                                          '/', chargedunit);
      if rankd_z2 = k then write (outfile, #09, controlline [81]);
      if rankn_z2 = k then write (outfile, #09, controlline [82],
                                          '/', neutralunit);
      if rankbalance2 = k then write (outfile, #09, controlline [83]);
      if rank_q = k then write (outfile, #09, controlline [84]);
      if rank_nq = k then write (outfile, #09, controlline [85],
                                        '/', backchargeunit);
      if rank_nd = k then write (outfile, #09, controlline [86],
                                        '/', backgroundunit);
      if ranktable_p = k then write (outfile, #09, header_p);
      if ranktable_n = k then write (outfile, #09, header_n);
      if ranktable_z = k then write (outfile, #09, header_z);
   end;
   writeln (outfile);
   for interval := 0 to numberofintervals do begin
      minute := interval * minutesininterval;
      moment := minute * momentsinminute;
      x := minute / timeunit;
      if timeunit < 15 then write (outfile, x:1:0) else
      if timeunit < 90 then write (outfile, x:4:2) else
                            write (outfile, x:5:3);
      effectivedialength := (backgrounddiameter [interval] - 1.5)
                           * backgroundconcentration [interval];
      for k := 1 to 99 do begin
         if ranknucl_p = k then writeresult2 (topnucleation_p *
                                              nucleationratio [moment]);
         if ranknucl_n = k then writeresult2 (topnucleation_n *
                                              nucleationratio [moment]);
         if ranknucl_z = k then writeresult2 (topnucleation_z *
                                              nucleationratio [moment]);
         if f_ranknucl_p = k then writeresult2 (f_topnucleation_p *
                                                nucleationratio [moment]);
         if f_ranknucl_n = k then writeresult2 (f_topnucleation_n *
                                                nucleationratio [moment]);
         if f_ranknucl_z = k then writeresult2 (f_topnucleation_z *
                                                nucleationratio [moment]);
         if rankion_p = k then writeresult2 (n_ion_p [interval]/ionunit);
         if rankion_n = k then writeresult2 (n_ion_n [interval]/ionunit);
         if rankd_p = k then writeresult2 (dp1 [interval]);
         if rankn_p = k then writeresult2 (np1 [interval]/chargedunit);
         if rankd_n = k then writeresult2 (dn1 [interval]);
         if rankn_n = k then writeresult2 (nn1 [interval]/chargedunit);
         if rankd_z = k then writeresult2 (dz1 [interval]);
         if rankn_z = k then writeresult2 (nz1 [interval]/neutralunit);
         if rankbalance1 = k then writeresult2 (balance1 [interval]);
         if rankd_p2 = k then writeresult2 (dp2 [interval]);
         if rankn_p2 = k then writeresult2 (np2 [interval]/chargedunit);
         if rankd_n2 = k then writeresult2 (dn2 [interval]);
         if rankn_n2 = k then writeresult2 (nn2 [interval]/chargedunit);
         if rankd_z2 = k then writeresult2 (dz2 [interval]);
         if rankn_z2 = k then writeresult2 (nz2 [interval]/neutralunit);
         if rankbalance2 = k then writeresult2 (balance2 [interval]);
         if rank_q = k then writeresult2 (100 * backcharge [interval]/
                                          backgroundconcentration [interval]);
         if rank_nq = k then writeresult2 (backcharge [interval]/
                                           backchargeunit);
         if rank_nd = k then writeresult2 (0.001*effectivedialength/
                                                 backgroundunit);
         if ranktable_p = k then begin
             for i := 1 to numberoffractions do
                            buffer [i] := table_p [i, interval];
             for i := 1 to numberoffractions do writeadjust (i);
         end;
         if ranktable_n = k then begin
             for i := 1 to numberoffractions do
                            buffer [i] := table_n [i, interval];
             for i := 1 to numberoffractions do writeadjust (i);
         end;
         if ranktable_z = k then begin
             for i := 1 to numberoffractions do
                            buffer [i] := table_z [i, interval];
             for i := 1 to numberoffractions do writeadjust (i);
         end;
      end; {of k}
      writeln (outfile);
   end; {of interval counter}
   closefile (outfile);
end; // of standard table

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{COMPILE and WRITE PARAMETER LIST}
if deliverparameterlist then begin
   s := nameofresult; s [1] := 'L'; s := s + '.txt';
   assignfile (outfile, s); rewrite (outfile);
   if delimiter = 0 then s := ' ' else
   if delimiter = 2 then s := ','
                    else s := #09;
   for i := 2 to 57 do write (outfile,
      controlline [i], '=', controlvalue [i], s);
   write (outfile,
      'sct=', numberofsections, s,
      'max=', finalsize:4:2, s,
      '?=', indeterminacy, s,
      'unit=', timeunit:1:0);
   if n_growthtable > 1 then write (outfile, ' tabulated growth');
   writeln (outfile);
   closefile (outfile);
   writeln;
end; // of parameterlist

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

{WRITE TABLE OF THE LAST DISTRIBUTION}
if deliverdistribution then begin
   s := nameofresult; s [1] := 'D'; s := s + '.xl';
   assignfile (outfile, s); rewrite (outfile);
   if diameterormobility = 1
      then writeln (outfile, 'd'#9'N+'#9'N-'#9'N0')
      else writeln (outfile, 'Z'#9'N+'#9'N-'#9'N0'#9'd'#9'Gr_fr');
   k := numberofintervals;
   for fraction := 1 to numberoffractions do begin
      if diameterormobility = 1 then i := fraction
                                else i := numberoffractions - fraction + 1;
      if diameterormobility = 1 then writeln (outfile, d_ctr [i]:0:3, #9,
         table_p [i, k]:0:2, #9, table_n [i, k]:0:2, #9, table_z [i, k]:0:1)
                                else writeln (outfile, z_ctr [i]:0:4, #9,
         table_p [i, k]:0:2, #9, table_n [i, k]:0:2, #9, table_z [i, k]:0:1, #9,
          d_ctr [i]:0:3, #9, 3600 * freeairgrowthrate (d_ctr [i], 1, 0):0:3);
   end;
   closefile (outfile);
   writeln;
end;

{COMPILE PLOTTABLE}
if deliverplottable then begin
  {make fraction diameter upper borders d_fr and centers}
   a := 0.5; b := 10;
   for fraction := 0 to 20 do
      d_fr [fraction] := exp (ln (a) + fraction * ln (b / a) / 20);
   for fraction := 1 to 20 do
      d_ctr [fraction] := sqrt (d_fr [fraction] * d_fr [fraction - 1]);

  {make distribution subtables for fraction concentrations

    pseudoindex:       x                     y
                --+-------+-------+-------+-------+---
    section index:        i      i+1     j-1      j    }

   for interval := 0 to numberofintervals do begin

      {distribute positive ions into buffer}
      a := 2 * sqr (ln (g_sigma_p)); b := 0;
      for fraction := 1 to 20 do begin
         d := sqr (ln (d_ion_p) - ln (d_ctr [fraction])) / a;
         if d > 22 then buffer [fraction] := 0
         else buffer [fraction] := exp (-d);
         b := b + buffer [fraction];
      end;
      for fraction := 1 to 20 do buffer [fraction] :=
         n_ion_p [interval] * buffer [fraction] / b;
      {add particles}
      for fraction := 1 to 20 do begin {positive particles}
         x := (d_fr [fraction - 1] - birthsize) / dd; {pseudoindex}
         y := (d_fr [fraction] - birthsize) / dd; {pseudoindex}
         i := trunc (x) + 1; // from
         j := trunc (y) + 1; // until
         if (i > maxlast) or (j < first) then c := 0 else begin
            if i < first then begin i := first; x := first end;
            if j > maxlast then begin j := maxlast; y := maxlast end;
            c := (i - x) * n_particle_p [i, interval]
               + (y - j + 1) * n_particle_p [j, interval];
            if (j - 1) > (i + 1) then
            for k := i + 1 to j - 1 do
               c := c + n_particle_p [k, interval];
         end;
         table_p [fraction, interval] := c + buffer [fraction];
      end {of positive particles};

     {distribute negative ions into buffer}
      a := 2 * sqr (ln (g_sigma_n)); b := 0;
      for fraction := 1 to 20 do begin
         d := sqr (ln (d_ion_n) - ln (d_ctr [fraction])) / a;
         if d > 22 then buffer [fraction] := 0
         else buffer [fraction] := exp (-d);
         b := b + buffer [fraction];
      end;
      for fraction := 1 to 20 do buffer [fraction] :=
         n_ion_n [interval] * buffer [fraction] / b;
      {add particles}
      for fraction := 1 to 20 do begin {negative particles}
         x := (d_fr [fraction - 1] - birthsize) / dd; {pseudoindex}
         y := (d_fr [fraction] - birthsize) / dd; {pseudoindex}
         i := trunc (x) + 1;
         j := trunc (y) + 1;
         if (i > maxlast) or (j < first) then c := 0 else begin
            if i < first then begin i := first; x := first end;
            if j > maxlast then begin j := maxlast; y := maxlast end;
            c := (i - x) * n_particle_n [i, interval]
               + (y - j + 1) * n_particle_n [j, interval];
            if (j - 1) > (i + 1) then
            for k := i + 1 to j - 1 do
               c := c + n_particle_n [k, interval];
         end;
         table_n [fraction, interval] := c + buffer [fraction];
      end {of negative particles};

      for fraction := 1 to 20 do begin {neutral particles}
         x := (d_fr [fraction - 1] - birthsize) / dd; {pseudoindex}
         y := (d_fr [fraction] - birthsize) / dd; {pseudoindex}
         i := trunc (x) + 1;
         j := trunc (y) + 1;
         if (i > maxlast) or (j < first) then c := 0 else begin
            if i < first then begin i := first; x := first end;
            if j > maxlast then begin j := maxlast; y := maxlast end;
            c := (i - x) * n_particle_z [i, interval]
               + (y - j + 1) * n_particle_z [j, interval];
            if (j - 1) > (i + 1) then
            for k := i + 1 to j - 1 do
               c := c + n_particle_z [k, interval];
         end;
         table_z [fraction, interval] := c;
      end {of neutral particles maybe used in future versions of the program};

   end; {of distribution subtables that include fraction concentrations}

  {WRITE PLOTTABLE}
   c := ln (d_fr [2] / d_fr [1]) / ln (10);
   s := nameofresult; s [1] := 'P'; s := s + '.xl';
   assignfile (outfile, s); rewrite (outfile);
   p := numberofintervals div 4;
   write (outfile, diatablecomposition);
   for i := 1 to 20 do write (outfile, #9, d_ctr [i]:0:3);
   for i := 1 to 20 do write (outfile, #9, d_ctr [i]:0:3);
   writeln (outfile);
   for k := 0 to p + numberofintervals do begin
      interval := k - p;
      write (outfile, interval);
      if interval < 0 then interval := 0;
      if diatablecomposition < 3 then for i := 1 to 20 do
                                 writeresult1 (table_p [i, interval] / c)
                                 else for i := 1 to 20 do
                                 writeresult1 (table_n [i, interval] / c);
      if diatablecomposition < 2 then for i := 1 to 20 do
                                 writeresult1 (table_n [i, interval] / c)
                                 else for i := 1 to 20 do
                                 writeresult1 (table_z [i, interval] / c);
      writeln (outfile);
   end;
   closefile (outfile);

end; {if deliverplottable}

   s := copy (nameofresult, 2, 3);
   write ('Files');
   if deliverstandardtable then write (' T', s);
   if deliverparameterlist then write (' L', s);
   if deliverdistribution then write (' D', s);
   if deliverplottable then write (' P', s);
   writeln (' are successfully saved.');
   writeln ('Processing time ', 86400 * (now - time0):1:2, ' seconds.');
   write ('Press ENTER for exit!'); readln;

except end;
end.

{Demo version of the burstcontrol.txt:

======= BURSTCONTROL_DEMO.TXT for BURSTSIMULATOR version HT20110306 =======
= The lines beginning with a symbol "=" in the first position and all lines
= after THE END are the comment lines. The lines beginning with a different
= symbol are the control lines. A control line consists of a control value
= followed by a space and a comment. A dual control value consists of the
= free-air value, slash, and the forest value. The number and the identifier
= of the control variable are shown at the end of a control line.
= The control values can be modified by the user of the simulator. 
=== NB: Do not add or delete any control lines in this file! ===
=== Symbol # indicates that the value must be a whole number.
=== Symbol * allows to present slash-separated dual values e.g. 4.5/6.7
=== where the first is valid in the free air and the second in the forest.
= Full instruction is available as a separate document.
===========================================================================
240 #min : Time period under consideration (<= 1440), 1
40 #: Number of evolution steps in a minute (recommended 20...60), 2:nst
1 : Acceleration coefficient dd/dd_min, normal value is 1, 3:acc
0 Celsius: Air temperature, 4:t
1013 millibar: Air pressure, 5:p
=== IONIZATION PARAMETERS ===
4/5 *cm-3 s-1: Initial ionization rate, 6:I0
4.5/5 *cm-3 s-1: Halftime ionization rate, 7:I1
5 *cm-3 s-1: Final ionization rate, 8:I2
1.4 cm2 V-1 s-1: Electric mobility of positive cluster ions, 9:Z+
1.6 cm2 V-1 s-1: Electric mobility of negative cluster ions, 10:Z-
1.5e-6 cm3 s-1: Cluster ion mutual recombination coefficient, 11:rec
2 g cm-3: Density of ions (standard value is 2.08), 12:idn
=== NUCLEATION PARAMETERS ===
1.5 nm: Birth size of particles, 13:d0
0 *cm-3 s-1: Maximum nucleation rate for positive ion nucleation, 14:J+
0.1 *cm-3 s-1: Maximum nucleation rate for negative ion nucleation, 15:J-
1/2 *cm-3 s-1: Maximum nucleation rate for neutral nucleation, 16:J0
10 min: Rise time of the nucleation activity, 17:rt
3 #: Shape code of rise: 1=linear, 2=sinus, 3=square_of_sinus, 18:shr
60 min: Time of steady nucleation activity, 19:st
30 min: Time of dropping nucleation activity, 20:dt
3 #: Shape of dropping: 1=linear, 2=sinus, 3=square_of_sinus, 21:shd
=== PARAMETERS OF THE FIRST CONDENSING SUBSTANCE ===
2 g cm-3: Density of growth units, 22:u1r
0.115 nm: Extra distance of the Van der Waals capture, 23:h
0.55 nm: Diameter of a growth unit, 24:u1d
149 ^3 : Polarizability, 25:pol1
0 : Nadykto-Yu dipole enhanchement factor for d = 1 nm, 26:Yu1
0 : Nadykto-Yu dipole enhanchement coefficient for d = 2 nm, 27:Yu2
1 *nm/h: Initial plain Knudsen growth rate of neutral particles, 28:G10
2/3 *nm/h: Halftime plain Knudsen growth rate of neutral particles, 29:G11
1 *nm/h: Final plain Knudsen growth rate of neutral particles, 30:G12
2.5 nm: Critical size of quantum rebound, 31:dq0
600 K: Extra temperature of quantum rebound, 32:Tq0
=== PARAMETERS OF THE SECOND CONDENSING SUBSTANCE ===
0.8 nm: Diameter of a growth unit, 33:u2d
32 ^3: Polarizability, 34:pol2
2 *nm/h: Initial plain Knudsen growth rate of neutral particles, 35:G20
3/7 *nm/h: Halftime plain Knudsen growth rate of neutral particles, 36:G21
3 *nm/h: Final plain Knudsen growth rate of neutral particles, 37:G22
3 nm: Critical diameter when the condensation starts, 38:nKd
2 : Power of nano-Koehler approximation, 39:nKp
=== PARAMETERS OF THE PRE-EXISTING BACKGROUND AEROSOL ===
50 nm: Initial average diameter of background aerosol particles, 40:bd0
60 nm: Halftime average diameter of background aerosol particles, 39:bd1
60 nm: Final average diameter of background aerosol particles, 42:bd2
1500 cm-3: Initial concentration of background aerosol particles, 43:bn0
2000 cm-3: Halftime concentration of background aerosol particles, 44:bn1
3000 cm-3: Final concentration of background aerosol particles, 45:bn2
=== PARAMETERS OF THE CONIFER FOREST ===
180 m: Air residence distance (time * wind) in forest, 46:rd
1 m s-1: Initial wind in the forest, 47:w0
1 m s-1: Halftime wind in the forest, 48:w1
1 m s-1: Final wind in the forest, 49:w2
0.9 mm: Conifer needle diameter, 50:dn
150 m-2: Conifer needle length in a unit volume m/m3, 51:L
=== PRESENTATION CONTROL ===
1 #: Particle distribution argument is 1) diameter or 2) mobility, 52:d/Z
1.333 nm or cm2 V-1 s-1: Low end of the particle presentation range, 53:dmn
7.5 nm or cm2 V-1 s-1: High end of the particle presentation range, 54:dmx
6 #: Number of size fractions in distribution tables (<= 99), 55:nfr
2 #: Scale of size or mobility 1) linear 2) logarithmic, 56:scl
3 #: Distr. 1) fraction concentrations 2) dN/dd 3) dN/d(log(d)), 57:dst
5 #min: Duration of a time interval in the output table, 58
1 #min: Unit of time in the output table, 59:unit
0 #: Delimiter in the signature file 0) gap 1) tab 2) comma, 60
TL : Output: T = st-table, L = par-list, D = distr, P = plot-table 61
=== RANKS OF VARIABLES AND SUBTABLES (RANK 0 = SKIP VARIABLE) ===
0 #: Free air positive nucleation rate, 62:J+a
0 #: Free air negative nucleation rate, 63:J-a
0 #: Free air neutral nucleation rate, 64:Joa
0 #: Forest positive nucleation rate, 65:J+f
1 #: Forest negative nucleation rate, 66:J-f
1 #: Forest neutral nucleation rate, 67:Jof
1 #: Concentration of positive ions, 68:n+
1 #: Concentration of negative ions, 69:n-
0 #: Average diameter of positive particles, 70:d+
1 #: Concentration of positive particles, 71:N+
0 #: Average diameter of negative particles, 72:d-
1 #: Concentration of negative particles, 73:N-
1 #: Average diameter of neutral particles, 74:do
0 #: Concentration of neutral particles, 75:No
0 #: Polarity balance index of charged particles, 76:balance5
0 #: Average diameter of positive particles between dmn&dmx, 77:d<+
0 #: Concentration of positive particles between dmn&dmx, 78:N<+
0 #: Average diameter of negative particles between dmn&dmx, 79:d<-
0 #: Concentration of negative particles between dmn&dmx, 80:N<-
0 #: Average diameter of neutral particles between dmn&dmx, 81:d<o
0 #: Concentration of neutral particles between dmn&dmx, 82:N<o
0 #: Polarity balance index between dmn&dmx, 83:balance<
0 #: Charge of background aerosol particles multiplied by 100, 84:100q
0 #: Charge concentration of background aerosol, 85:Nq
0 #: Effective diameter concentration of background aerosol, 86:Nd
1 #: Distribution subtable of positive particles, 87:p
1 #: Distribution subtable of negative particles, 88:n
0 #: Distribution subtable of neutral particles, 89:z
=== SCALE DENOMINATORS FOR DIAGRAMS ===
100 #: Concentration of ions cm-3, 90
50 #: Concentration of charged nanoparticles cm-3, 91
500 #: Concentration of neutral nanoparticles cm-3, 92
-10 #: Concentration of background aerosol charge m-3, 93
100 #: Effective diameter concentration of background aerosol m-2, 94
=== INFORMATION ABOUT PLOT TABLE ===
1.18 : geometric std of the distribution of positive ions, 95:gs+
1.22 : geometric std of the distribution of negative ions, 96:gs-
1 #: plot table composition: 1 = + & -, 2 = + & neutral, 3 = - & neutral, 97
=== REPLACEMENT FOR INDETERMINED VALUES ===
? : Symbol(s) of indeterminacy (- = gap, # = presumable limit), 98
=== THE END ===
COMMENTS:
1. The columns of the output standard table are delivered in the order
   of ranks and when of equal rank, then in the order of control lines.
   The columns of zero rank are skipped in the output table.
2. Scale denominators are applied to the integral concentrations only,
   the distributions are delivered in cm-3 or cm-3 nm-1 independent
   of denominator values.

