diff --git a/src/cstream.c b/src/cstream.c
new file mode 100644
index 0000000000000000000000000000000000000000..a91857f6f1f259caec51311a251ff219b8d0d030
--- /dev/null
+++ b/src/cstream.c
@@ -0,0 +1,459 @@
+/** @file cstream.c
+ *  @brief General initialization/termination functions.
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <time.h>
+#define ALLOC_PARAMS
+
+#include "cstream.h"
+#include "parameters.h"
+#include "proto.h"
+#include "species.h"
+
+double E_x, E_y, E_z;
+/**< For time-dependent electric fields E_x, E_y, E_z there is a difference between
+   E0 (the highest field) and E (the field at a given time). */
+
+extern double pois_inhom_fixed_q_t;
+/**< The same is true for time-dependent fixed charges (useful to model
+   the charging of a cloud in modeling sprites. */
+
+double *dr, *dz, dtheta;
+/**< dr[i] and dz[i] represent the grid size at level @a i.  Note that
+   @a i can go from -extra_pois_levels...max_levels. */
+double *dr_start, *dz_start;
+
+double *w2k, *wk;
+/**< w2k[k] = 2 / dtheta^2 (1 - cos(k dtheta)). This numbers appear due to the
+   finite-differences derivatives of a Fourier transform.  If we would make
+   a continuous fourier transform, they will approach k^2.
+*/
+
+
+char *invok_name = "cstream";
+/**< Name of this program. */
+
+const double twopi = 6.283185307179586477L;
+/**< \f$ 2 \times \pi \f$. */
+
+const double invfourpi = 0.079577471545947667884441882L;
+/**< \f$ 1 / (4 \pi )\f$ */
+
+const double invpi32 = 0.179587122125166561689081L;
+/**< \f$ 1 / (\pi ^ {3/2} )\f$ */
+
+/** @brief Initializes the parameters.
+ *
+ *  The global parameters will get their default values.
+ *  These values correspond with the values in input/defaults.cfg
+ */
+void
+init_parameters(void)
+{
+    debug (1, "entry init_parameters\n");
+    /* An identification name for this run */
+    prog_id = "example";
+
+    /* Output directory.  */
+    output_dir = "output/";
+
+    /* Kinetics input file */
+    kin_input = "input/kinetic_example.cfg";
+
+    /* If restart is TRUE, the simulation will continue with data from a previous run */
+    restart = 0;
+
+    /* If restart is TRUE, the name of the file with data from previous run, otherwise empty */
+    load_file = "";
+
+    /* Time interval for output to be written to disk */
+    output_dt = 0.100;
+
+    /* Output of the Poisson grids, including the potential? */
+    pois_output = 0;
+
+    /* Margin outside the grids in the output of the cdr equation */
+    cdr_output_margin =  0;
+
+    /* Margin outside the grids in the output of the poisson equation */
+    pois_output_margin = 1;
+
+    /* If the time steps are smaller than this number, the program issues a warning */
+    warn_min_timestep = 1e-06;
+
+    /* Maximum disk space, in Mb, to use */
+    max_disk_space_mb = 1048576;
+
+    /* Number of R and Z gridpoints at level 0 */
+    gridpoints_r = 600;
+    gridpoints_z = 600;
+
+    /* Number of azimuthal gridcells and modes */
+    max_ntheta = 1;
+
+    /* Initial and end time */
+    start_t = 0.0;
+    end_t = 0.12;
+
+    /* Attempted timestep.  The actual timestep may be larger */
+    attempt_dt = 50.0;
+
+    /* Extra levels for the Poisson solver */
+    extra_pois_levels = 2;
+
+    /* Maximum level of refinement. Use a big number here */
+    max_levels = 64;
+
+    /* Error threshold that leads to refinement in the Poisson code. */
+    pois_max_error = 0.001;
+
+    /* Maximum level of refinement in the Poisson equation. */
+    pois_max_level = 3;
+
+    /* Extra levels for the photo-ionization solver */
+    extra_photo_levels = -1;
+
+    /* Maximum level of refinement in the photo-ionization solver. */
+    photo_max_level = 4;
+
+    /* Error threshold that leads to refinement in the photo-ionization code. */
+    photo_max_error = 0.01;
+
+    /* Photo-ionization boundary condition at r = L_r, z = 0, z = L_z.  1 for Hom. Neumann, -1 for Hom. Dirichlet */
+    photo_bnd_right  = -1;
+    photo_bnd_bottom = -1;
+    photo_bnd_top    = -1;
+
+    /* Extra levels for the photo-ionization solver */
+    extra_photo_levels_2 = -1;
+
+    /* Maximum level of refinement in the photo-ionization solver. */
+    photo_max_level_2 = 4;
+
+    /* Error threshold that leads to refinement in the photo-ionization code. */
+    photo_max_error_2 = 0.01;
+
+    /* Photo-ionization boundary condition at r = L_r, z = 0, z = L_z.  1 for Hom. Neumann, -1 for Hom. Dirichlet */
+    photo_bnd_right_2  = -1;
+    photo_bnd_bottom_2 = -1;
+    photo_bnd_top_2    = -1;
+
+    /* Particles boundary condition at z = 0, z = L_z, r = L_r.  1 for Hom. Neumann, -1 for Hom. Dirichlet */
+    cdr_bnd_bottom = 1;
+    cdr_bnd_top    = 1;
+    cdr_bnd_right  = 1;
+
+    /* Potential boundary condition at r = L_r, z = 0,  z = L_z.  1 for Hom. Neumann, -1 for Hom. Dirichlet */
+    pois_bnd_right  = -1;
+    pois_bnd_bottom = -1;
+    pois_bnd_top    = -1;
+
+    /* Maximum advection and diffusion Courant number */
+    nu_a = 0.2;
+    nu_d = 0.2;
+
+    /* Maximum ratio of dt/relaxation time */
+    nu_rt = 0.2;
+
+    /* Maximum ratio of change of the densities (set to a very large number to ignore) */
+    nu_f = 1e+20;
+
+    /* Refinement threshold for the electric field */
+    ref_threshold_eabs = 0.2;
+
+    /* Maximum refinement level reached through ref_threshold_eabs */
+    ref_level_eabs = 4;
+
+    /* Refinement threshold for the curvature of the charge, densities */
+    ref_threshold_charge = 0.004;
+    ref_threshold_dens = 0.004;
+
+    /* Refinement threshold for the densities in the leading edge */
+    ref_threshold_edge = 10000.0;
+
+    /* r-length and  z-length  of the minimal refinement area in the cdr equation */
+    cdr_brick_dr = 8;
+    cdr_brick_dz = 8;
+
+    /* Maximum level of refinement in the Fluid equation. */
+    cdr_max_level = 3;
+
+    /* Interpolation method for the grid interior, and grid boundaries (0=zero_masses, 1=quadratic_masses [default], 2=wackers_masses, 3=quadlog */
+    cdr_interp_in = 1;
+    cdr_interp_bnd = 1;
+
+    /* Length in r and z of the complete domain */
+    L_r = 13044.0;
+    L_z = 13044.0;
+
+    /* Isotropic difussion coefficient */
+    diffusion_coeff = 0.1;
+
+    /* Whether the code includes photoionization or not */
+    has_photoionization = 1;
+
+    /* The name of a file from which we can read the photoionization parameters */
+    photoionization_file = "input/air760torr.photo";
+
+    /* Rate of dissociative attachment */
+    attachment_rate = 0.0;
+
+    /* E0 in the exp(-E0/E) factor in the attachment expression. */
+    attachment_E0 = 0.0;
+
+    /* x-, y- and z-component of the external electric field */
+    E0_x = 0.0;
+    E0_y = 0.0;
+    E0_z = -0.06;
+
+    /* Rise time of the electric field (0 for instantaneous rise) */
+    rise_time = 0.0;
+
+    /* Time to switch off the electric field (0.0 means never) */
+    off_time = 0.0;
+
+    /* x-, y- and z-width of the initial seed */
+    seed_sigma_x = 0.0;
+    seed_sigma_y = 0.0;
+    seed_sigma_z = 0.0;
+
+    /* Number of electrons in the initial seed */
+    seed_N = 0.0;
+
+    /* Initial at z=0 densities of electrons and ions */
+    background_ionization = 0.0;
+
+    /* Length of exponential increase of the pre-ionization (for atmospherical models) */
+    background_increase_length = 0.0;
+
+    /* Use the point-plane geometry? */
+    pois_inhom = 1;
+
+    /* Number of mirror charges to use */
+    pois_inhom_reflections = 4;
+
+    /* Length and radius of the needle */
+    needle_length = 2500.0;
+    needle_radius = 400.0;
+
+    /* If nonzero, the charge is fixed, not floating (simulation of charged clouds close to the earth surface) */
+    pois_inhom_fixed_q = 0.0;
+
+    /* Constant ionization rate */
+    constant_source = 0.0;
+
+    /* Initial perturbation to the axisymmetric configuration */
+    perturb_epsilon = 0.0;
+
+    /* Perturb only modes up to perturb_max_k (large number to perturb all) */
+    perturb_max_k = 1024;
+
+    /* 1 if the sprite module is activated, 0 otherwise */
+    sprite_module = 0;
+
+    /* Lenght of exponential decay of the density w/r to altitude */
+    dens_decay_len = 0.0;
+
+    /* Density at z = 0 */
+    sprite_dens_0 = 0.0;
+
+    /* Quenching density */
+    sprite_dens_q = 0.0;
+
+    /* Sign of the sprite head that we are following (the other will not be reliable */
+    sprite_sign = -1;
+    debug (1, "exit init_parameters\n");
+}
+
+  int n, i;
+/** @brief Initializes the grid sizes.  The parameters dr_root and dz_root
+ *  specify the grid size at level 0 (dz_root = dz[0], dr_root[0] = dr[0])
+ */
+static void
+init_gridsizes_a (void)
+{
+  int n, i;
+  double root_dr, root_dz;
+
+  root_dr = L_r / gridpoints_r;
+  root_dz = L_z / gridpoints_z;
+
+  n = max_levels + extra_pois_levels + 1;
+
+  dr_start = (double *) xmalloc(sizeof(double) * n);
+  dz_start = (double *) xmalloc(sizeof(double) * n);
+
+  dr = dr_start + extra_pois_levels;
+  dz = dz_start + extra_pois_levels;
+
+  dr_start[0] = root_dr * (1 << extra_pois_levels);
+  dz_start[0] = root_dz * (1 << extra_pois_levels);
+
+  for (i = 1; i < n; i++) {
+    dr_start[i] = dr_start[i - 1] / 2.0;
+    dz_start[i] = dz_start[i - 1] / 2.0;
+
+    debug (3, "dr[%d] = %e\n", i - extra_pois_levels,
+	   dr[i - extra_pois_levels]);
+    debug (3, "dz[%d] = %e\n", i - extra_pois_levels,
+	   dz[i - extra_pois_levels]);
+  }
+
+  dtheta = twopi / max_ntheta;
+}
+
+/** @brief Initializes the vector of wk's for the Helmholtz equation. */
+static void
+init_wk_a (void)
+{
+  int k;
+  double twobydtheta2;
+
+  debug (2, "init_w2k_a()\n");
+
+  w2k = (double*) xmalloc (sizeof(double) * max_ntheta);
+  wk  = (double*) xmalloc (sizeof(double) * max_ntheta);
+
+  twobydtheta2 = 2. / (dtheta * dtheta);
+
+  for (k = 0; k < max_ntheta / 2 + (max_ntheta % 2); k++) {
+    double w2k_ = twobydtheta2 * (1 - cos (k * dtheta));
+    double re_wk, im_wk;
+
+    re_wk = (cos (k * dtheta) - 1) / dtheta;
+    im_wk = sin (k * dtheta) / dtheta;
+
+    wk[k] = re_wk;
+    w2k[k]= w2k_;
+
+    if (k != 0) {
+      wk[max_ntheta - k] = im_wk;
+      w2k[max_ntheta - k] = w2k_;
+    }
+  }
+
+  if ((max_ntheta % 2) == 0) {
+    w2k[max_ntheta / 2] = twobydtheta2 * (1 - cos (k * dtheta));
+    wk[max_ntheta / 2] = (cos (k * dtheta) - 1) / dtheta;
+  }
+
+  assert (wk[0] == 0.);
+}
+
+/** @brief Frees the space allocated for dr and dz. */
+static void
+free_gridsizes (void)
+{
+  free (dr_start);
+  free (dz_start);
+}
+
+/** @brief Frees the allocated space for \f$ \lambda \f$. */
+static void
+free_wk (void)
+{
+  free (wk);
+  free (w2k);
+}
+
+/** @brief  Here all the initialization calls. */
+void
+cstream_init (void)
+{
+  /* For the perturbations, it is better to start with different seeds at
+     each run. */
+  srand (time (0));
+
+  init_gridsizes_a ();
+  init_wk_a ();
+  react_init ();
+  cdr_init ();
+  pois_init ();
+  photo_init ();
+  if (has_photoionization) {
+    printf("Photoionization\n");
+    photo_load_file (photoionization_file);
+  } else {
+    printf("No photoionization\n");
+  }
+
+  if (sprite_module) spr_init ();
+}
+
+/** @brief Here all the cleaning and finishing calls. */
+void
+cstream_end (void)
+{
+  photo_unregister_all ();
+  free_gridsizes ();
+  free_wk ();
+  cdr_end ();
+}
+
+
+/** @brief When we use a time-varying external field, we update the field
+ *  components at each timesteps. */
+void
+cstream_set_field_at_time (double t)
+{
+  static int rise_reached = FALSE, off_reached = FALSE;
+  double factor;
+
+  if (rise_reached && off_reached)
+    return;
+
+  if (t >= rise_time) {
+    rise_reached = TRUE;
+  }
+
+  if (t >= off_time || off_time == 0.0) {
+    off_reached = TRUE;
+  }
+
+  /* We make sure that we never use a field _larger_ than E0 and that
+     if rise_time is 0 we do not calculate t / 0.0.  */
+  factor = t >= rise_time? 1.0: (t / rise_time);
+
+  factor = (off_time > 0.0 && t >= off_time)? 0.0: factor;
+
+  E_x = factor * E0_x;
+  E_y = factor * E0_y;
+  E_z = factor * E0_z;
+
+  if (pois_inhom_fixed_q != 0.0) {
+    pois_inhom_fixed_q_t = factor * pois_inhom_fixed_q;
+  }
+
+  return;
+}
+
+/** @brief Functions for a constant external electric field given in (x, y, z) components. */
+double
+e0_r (double r, double z, double theta)
+{
+  return E_x * cos (theta) + E_y * sin (theta);
+}
+
+/** @brief Functions for a constant external electric field given in (x, y, z) components. */
+double
+e0_z (double r, double z, double theta)
+{
+  return E_z;
+}
+
+/** @brief Functions for a constant external electric field given in (x, y, z) components. */
+double
+e0_theta (double r, double z, double theta)
+{
+  return -E_x * cos (theta) + E_y * sin (theta);
+}
+
+/** Initializes the component of the external electric field in \f$r\f$-direction. */
+decl_field_comp(r) = e0_r;
+/** Initializes the component of the external electric field in \f$z\f$-direction. */
+decl_field_comp(z) = e0_z;
+/** Initializes the component of the external electric field in \f$\theta\f$-direction. */
+decl_field_comp(theta) = e0_theta;