The EIRENE Code User Manual

including: B2-EIRENE interface

Version: 11/2009

D. Reiter

EIRENE

THE GREEK GODDESS OF PEACE

The Greek author Hesiodos wrote a genealogy of the gods, The Theogony, in the 8th century B.C.. According to him, EIRENE and her sisters, Eunomia and Dike, were the daughters of Zeus and Themis. These sisters were called the horae, the Greek word for the right time, hora.

Abstract

The EIRENE neutral gas transport code is described. This code resorts to a combinatorial discretization of general 3 dimensional computational domains. It is a multi-species code solving simultaneously a system of time dependent (optional) or stationary (default) linear kinetic transport equations of almost arbitrary complexity. A crude model for transport of ionized particles along magnetic field lines is also included. EIRENE is coupled to external databases for atomic and molecular data and for surface reflection data, and it calls various user supplied routines, e.g. for exchange of data with other (fluid-) transport codes. The main goal of code development was to provide a tool to investigate neutral gas transport in magnetically confined plasmas. But, due to its flexibility, it also can be used to solve more general linear kinetic transport equations, by applying a stochastic rather than a numerical or analytical method of solution. In particular, options are retained to reduce the model equations to the theoretically important case of the one speed transport problem.
Major applications of EIRENE are in connection with plasma fluid codes, in particular with the various versions of the B2 code [1]. The semi-implicit iterative coupling method of B2-EIRENE [2], [3] and it's implementation (code segment: EIRCOP) are also described.

Foreword to 2nd edition

The first edition of this EIRENE code user manual was published as KFA report JÜL-2599 in March 1992, [4]. It basically was a collection of my notes on the meaning of the various input data options for EIRENE. Over the years too many such input flags had accumulated to memorize the meaning of each individual one.
Since the EIRENE code became a quite popular tool also for many other users, it was decided to provide these notes as a kind of user's manual in a somewhat completed and edited form.
In the meantime the distribution of the EIRENE code has become even wider, and it seems timely to update the previous manual, although most of it's content is still a relevant source of information for a user of the code. Apart from several minor corrections, e.g., of spelling errors and unclear language, the major new features as compared to the previous edition are the following:

- time dependent mode + snapshot estimators (section 2.13)
- internally consistent Ï-integral" approach for coupled neutral (kinetic) - plasma (fluid) simulations
- individual background ïon" temperatures TI(IPLS) for each background species (e.g., for neutral background particles)
- simulation of self collisions in BGK approximation
- some new collision processes and data, such as net re-combination energy losses (radiation + potential), elastic neutral - ion collisions, multi-step brake-up of molecules
- RAPS graphics interfaces
- user defined geometry block. In this new mode of operation, LEVGEO=6, EIRENE knows nothing (not even the dimensionality) about the geometrical aspects of the problem. Everything concerning grids, cell volumes, flight times within cell etc. is in user supplied routines. No default graphics options are available, of course, in this mode of operation.
- Particle tracing routines FOLATM and FOLMOL for neutral atoms and neutral molecules are combined into one routine FOLNEUT for neutral particles.
The particle tracing routine FOLION for charged test particles has been updated considerably, including now a simple approximation to the Fokker Planck collision operator (a "Krook collision operator").

An increasing number of EIRENE users is running cases in which EIRENE is coupled to a fluid plasma model to provide sources and sinks. The widely used B2-EIRENE code system [2] is one such example. A new section describing this part of the EIRENE code has been added. It describes the "sandwich" file EIRCOP which was written to permit linkage between EIRENE and a plasma fluid model (one, two, or three dimensional. Initial and boundary value problems can be treated in a rather self-consistent manner.
This code segment was developed largely under support of a KFA-EURATOM contract ([5]), and first applications have been published for ITER configurations, see [2] and [6]. Subsequently a large number of application runs, in particular those carried out at IPP-Garching (Ralf Schneider et al.) and AEA-Culham Laboratories (Geoff Maddison) have lead to identification of many critical issues and limitations of the code and to permanent improvements until these days. Typical convergence behaviour of the combined code system ("saturated residuals") is described in the new section 1.8, written in collaboration with G.P.Maddison, AEA Culham Laboratories.
Detlev Reiter, Spring 1998

Foreword to 3rd edition

During the year 2001 a major revision of the EIRENE code has been carried out, largely in order to implement a somewhat more modern FORTRAN structure into the code. EIRENE code versions since then have sometimes been referred to as "EIRENE-FACELIFT".
In particular a dynamical allocation of storage (rather than the hitherto necessary pre-assignment of storage in the PARAMETER statements collected in the previous PARMUSR file) is now in place. This dynamic memory managment has led to a replacement of the previous "Common-Blocks" now by "Modules", and the entire elimination of the "Equivalence" statements used for storage economy in earlier versions.
Further main upgrades are related to:
• incident-species dependent surface interaction models, pumping speeds, sputter models etc., can now be controlled by input flags, rather than the previous quite difficult implementation via the USR-routines.
• integration volumes of volume averaged tallies can now be larger than the grid cell volumes, i.e., the plasma background and geometrical discretisation, which is not very storage sensitive, can be made much finer than the neutral particle volume averaged tally profile output.
• discretisation of general multiply connected 3D volumes by tetrahedron-grids.
• photons as new type of test particle species. The photon_dummy module allows simple photon transport (one-speed Boltzmann-equation) with spatial profiles of spontaneous emission, Doppler line broadening, for optically thin lines only, with inclusion of (multiple) surface reflection. The full photon module additionally allows for various other line broadening mechanisms and photon re-absorption, stimulated emission, radiation trapping (excitation-hopping), iteration with neutral gas excitation population kinetics, etc...
EIRENE versions 2003 and younger have been compiled and successfully tested at FZ-Jülich on the following systems and compilers:
Linux: Suse 11.1 and all predecessors down to Suse 6.***:
• pgf90 Portland Compiler Version 5.2-4 - 10.1
• ifort Intel Fortran Compiler Version 8.1 - 11.1
• lf95 Lahey Fortran Compiler Version 6.2
• nagf95 NAG Fortran Compiler Version 5.0 - 5.2
AIX: AIX 4.3 und 5.2
• xlf95 XL Fortran for AIX Version 8.1
Windows: Windows 2000 und Windows XP
• Compaq Visual Fortran Version 6.6
Detlev Reiter, Spring 2010

Contents

1  The neutral gas transport equation;
Monte-Carlo terminology

1.1  The linear Boltzmann equation for the distribution function f
1.2  The linear integral equation for the collision density Ψ
1.2.1  The Green function concept
1.3  Monte Carlo solution of equation 1
1.3.1  Unbiased estimators
1.3.2  Scaling of tallies
1.3.3  Efficiency (FOM)
1.3.4  Source sampling
1.3.5  Sampling from the transport kernel T
1.3.6  Sampling from the collision kernel C
1.3.7  elastic collisions
1.3.8  charge exchange
1.3.9  electron-impact collisions
1.3.10  general heavy-particle-impact collisions
1.3.11  photon processes (emission, absorption, scattering
1.4  Surface Reflection Models
1.4.1  The Behrisch matrix reflection model
1.4.2  TRIM code database reflection models
1.5  Surface sources
1.5.1  The electrostatic sheath
1.6  Combinatorial description of geometry
1.7  time dependent mode of operation
1.8  non-linear effects: coupling to plasma fluid models
1.9  non-linear effects: neutral-neutral collisions
1.9.1  BGK approximation
1.9.2  Direct Simulation (DMCS) of self-collisions
1.10  Radiation transport, photon gas simulations
1.10.1  Line shape options
1.11  Charged Particle Transport
1.11.1  Orbit integration
1.11.2  Coulomb Collision Models
1.11.3  The full Fokker-Planck Collision Model, A
1.11.4  The full Fokker-Planck Collision Model, B: binary Monte Carlo collision procedure
1.12  EIRENE flow charts
2  Description of formatted input file
2.1  Input data for operating mode
2.1.1  The NLERG option for cell volumes
2.1.2  The NLMOVIE option for making movies of trajectories
2.2  Input for Standard Mesh
2.2.1  Mesh Parameters
2.3  The Input Block for Surfaces
2.3.1  The Input Block for "Non-default Standard Surfaces"
2.3.2  Input Data for "Additional Surfaces"
2.4  Input Data for Species Specification and Atomic Physics Module
2.4.1  Default atomic and molecular data
2.4.2  Neutral-Neutral collisions in BGK approximation
2.4.3  Fitting expressions (IFTFLG)
2.5  Input for Plasma Background
2.5.1  Derived Background Data
2.6  Input Data for Surface Interaction Models
2.7  Input data for Initial Distribution of Test Particles
2.7.1  Piecewise constant "Step-functions" for sampling
2.8  Additional Data for some Specific Zones
2.9  Data for Statistics and non-analog Methods
2.10  Data for additional volume and surface averaged tallies
2.10.1  Additional volume averaged tallies, tracklength estimator
2.10.2  Additional volume averaged tallies, collision estimator
2.10.3  Algebraic expression in volume averaged tallies
2.10.4  Algebraic expression in surface averaged tallies
2.10.5  Algebraic expression in surface averaged tallies
2.10.6  Energy spectra in selected cells or surfaces
2.11  Data for numerical and graphical Output
2.12  Data for Diagnostic Module
2.12.1  Line of sight: charge exchange spectrum
2.12.2  Line of sight: line emissivity
2.12.3  Line of sight: line shape
2.12.4  Line of sight: user defined integral
2.13  Data for nonlinear and time dependent Options
2.14  Data for interfacing Subroutine "INFCOP" (example)
2.14.1  Version B2-EIRENE-1999 and older
2.14.2  Version B2-EIRENE-2000 and younger
3  Problem specific Routines
3.1  Parameter Statements (for EIRENE-2001 or older)
3.2  The "Additional Tally" routines UPTUSR, UPCUSR, UPSUSR and UPNUSR
3.2.1  Tracklength estimated volume tallies, UPTUSR
3.2.2  Collision estimated volume tallies, UPCUSR
3.2.3  Snapshot estimated volume averaged tallies, UPNUSR
3.2.4  Surface averaged tallies, UPSUSR
3.3  The user surface reflection model REFUSR
3.4  The user source sampling routine SAMUSR
3.5  The user routines to overrule input data
3.5.1  The user geometry data routine GEOUSR
3.5.2  User supplied background data routine PLAUSR
3.6  The user routines for profiles PROUSR
3.7  User supplied post-processed tally routine TALUSR
3.8  User supplied "general geometry block"
3.8.1  Subroutine INIUSR
3.8.2  Subroutine LEAUSR
3.8.3  Subroutine TIMUSR
3.8.4  Subroutine VOLUSR
3.8.5  Subroutine NORUSR
4  Routines for interfacing with other codes:
EIRCOP

4.1  Routine for interfacing INFCOP
4.1.1  entry IF0COP
4.1.2  entry IF1COP
4.1.3  entry IF2COP(ISTRA)
4.1.4  entry IF3COP(ISTRAA,ISTRAE)
4.1.5  entry IF4COP
4.2  Routines for cycling of EIRENE with external codes: EIRSRT
4.3  Routines for special tallies needed for code coupling: UPTCOP
4.4  Statistical noise in Monte Carlo terms for external code, noise-residuals: STATIS_COP
5  Default EIRENE tallies, and selected Modules
5.1  Tables of EIRENE tallies
5.1.1  Current status, incl. photon gas tallies (Eirene-02 and younger)
5.1.2  old version, w/o photon gas tallies (Eirene-01 and older)
References
Introduction and General Information

This manual describes the input required by the EIRENE code to run a Monte Carlo study for a fully 3 dimensional simulation of linear transport (i.e., of test particles) in a prescribed background medium. Although the geometry of the problem and the interaction between test particle species and the background are in principle not subject to any restrictions, the aim of code development was to provide a tool for investigating neutral gas transport in tokamak plasmas. Consequently the choice of preprogrammed options has been made mainly with this application in mind. A large variety of problems in this field can be run without having to resort to any problem specific routines but instead by an appropriate setting of logical and numerical input flags.
Any user of EIRENE should be aware that this code is a moving target as is this manual.
Therefore it is possible that there are some inconsistencies between this description and a particular version of the code. The user should always first check subroutine INPUT, where most of the data are read in using hard wired formats, as most input errors will lead to a rapid exit in the initialization phase of a run.
This manual was written by the author of the code who often may not have been able to anticipate difficulties in understanding the use of EIRENE. He therefore gratefully acknowledges any suggestions to make this manual more informative and clear than it might be at the present status.
This code description consists of the following parts:
• In the first part an introduction is given to the general linear transport problem and its solution by Monte Carlo methods. Most of the terminology used in later sections is introduced there.
• In part two a description of the formatted input file required by EIRENE to run on a specific problem is given. It mainly consists of explanations of the meanings of the various input flags.
• In part three the most important problem specific routines are explained. At present we have restricted this part to routines for evaluation of "user requested tallies", namely the subroutines UPTUSR, UPCUSR and UPSUSR. Other often applied routines such as SAMUSR (user supplied source sampling distribution) or REFUSR (user supplied surface interaction model) are briefly described.
• In part four the package EIRCOP for interfacing EIRENE with other codes (e.g., the B2-EIRENE package) is described. Here mainly the location of the storage on the EIRENE work array RWK for plasma data and geometrical information is given. Also the implementation of the method of (semi-implicit) corrections (see section 1.8) at each plasma code time-step to the terms transferred from EIRENE to plasma fluid transport codes is described here.

Chapter 1 The neutral gas transport equation; Monte-Carlo terminology

General Remarks

To introduce the terminology used throughout this report, we briefly recall the basic definitions and principles of a Monte Carlo linear transport model, following the lead of many textbooks on Monte Carlo methods for computing neutron transport (see e.g., [7]). We begin with the linear transport equation for the pre-collision density, written as integral equation (linear non-homogeneous Fredholm integral equation of 2nd kind). Distinct from standard terminology in the (analytic) transport theory we do not discuss analytic properties of the various terms in the equation, but, instead, point out their probabilistic interpretation, as needed for a Monte-Carlo solution of that equation. Next (section 1.3) we sketch the Monte Carlo procedure to solve such equations, by referring to the two most often applied techniques: "track-length - and collision based estimators". In the third subsection we briefly describe the treatment of boundary conditions (models for interaction of the particles with surrounding surfaces) and discuss some special models which are in use for the neutral gas transport in fusion plasma devices. Then (section 1.5) we discuss the most important source function (non-homogeneous part of the integral equation) and its implementation in a Monte Carlo algorithm, namely the surface source of neutral particles due to recombination of ions incident on solid surfaces at the boundary or inside of the computational area. In section 1.6 we comment on the implementation of geometry within the framework of a Monte Carlo code in general terms and for the EIRENE code in particular. In section 1.7 the time dependent mode of operation of the EIRENE code is described. It merely amounts just to an increase in the dimensionality by one, by adding a time co-ordinate and treating it, formally, in a rather symmetric fashion with the other spatial coordinates. See also reference [8].
Two kinds of non-linearities may be accounted for:
In section 1.8 the nonlinear behavior resulting from background data, which depend on neutral particle transport (sources and sinks), is discussed. The algorithm of the B2-EIRENE code system [3] is described. More details on this can be found in the report [9].
In the final section of this introductory chapter 1.9 the non-linear BGK formalism and the direct simulation method (DMCS) for self collisions between neutral particles, as implemented in the EIRENE code, is described.

1.1  The linear Boltzmann equation for the distribution function f

The term μ-space in this report refers to the phase space of a single particle. The quantity of interest is then the one particle distribution function f
 f(r, v, i, t)     or     f(r, E, Ω, i, t)    or     f(x) ,
where the state x of the μ-space is characterized by a position vector r, a velocity vector v, a species index i (i stands for, e.g., H,D,T,D2,DT,He, CHn,...) and the time t.
The number density ni(r) for species i then reads:
 ni(r,t) = ⌠⌡ d3v  f(r,v,i,t)

Instead of v we sometimes utilize the kinetic energy E, E=m/2 ·v2, and the unit (speed) vector Ω = v/|v| in the direction of particle motion. Hence:
 f(r, v, i ,t) = ⎛⎝ m v ⎞⎠ ~ f (r, E, Ω, i, t)

We start by assuming that collisions are discontinuous events.
Further: Lets consider only one specific species i0, now omitting this species index. We assume that there are only collisions of this species i0 with only one further species (labelled b), and that exactly one particle of each of these species will also be present after the collision event (i.e., chemical reactions excluded for the time being, but see below). The familiar Boltzmann equation for the distribution function f for this species i0 reads

 [ ∂ ∂t + v·∇r + F(r,v,t) m ·∇v] f(r,v,t)
 =
 ⌠⌡ ⌠⌡ ⌠⌡ σ(v′,V′;v,V)|v′−V′| f(v′)fb(V′)
 −
 ⌠⌡ ⌠⌡ ⌠⌡ σ(v,V;v′,V′)|v−V| f(v) fb(V)
 +
 Q(r,v,t)
(1.1a)
Q is any external source (particles per unit time injected per unit volume in phase space).   Integrations are over the velocities v′,V and V′. σ(v′,V′;v,V) is the cross section for a binary particle collision process. It is equal to zero unless the conservation laws for total energy and momentum are fulfilled.
The first two arguments in σ, namely the velocities v′,V′ in the first integral, correspond to the species i0 and b, respectively, prior to a collision. These are turned into the post collision velocities v, V, again for species i0 and b, respectively. The first integral, therefore, describes transitions (v′,V′→ v,V) into the velocity space interval [v,v+dv] for species i0, and the second integral describes loss from that interval for this species.
Furthermore, m is the particle mass and F(r,v,t) is the volume force field.
The right hand side is the collision integral [(δf)/(δt)]|b. If there are more than just one possible type of collision partners, then the collision integral has to be replaced by a sum of collision integrals over all collision partners b, including, possibly, b=i0 (self collisions)
 δf δt |coll = δf δt |gain− δf δt |loss = ∑ b δf δt |collb.
(1.2)
This is readily generalized to the semi-classical Boltzmann equation for chemical reactions (including, for example, vibrational relaxation or exchange of internal energy as special cases) symbolized as i0 + j0 \rightleftarrows i1 + j1. These species indices label both the chemical species and/or the internal quantum state. In this case the sum in the collision integral is over j0, i1 and j1 and the cross sections in the corresponding collision integrals σi0j0i1j1(v,V,v′,V′) are differential for scattering at a certain solid angle and post collision energies with simultaneous transition from (i0, j0) to (i1, j1). Further generalizations to include particle splitting, absorption or fragmentation into more than 2 post collision products are straight forward, but can more conveniently be formulated in the C-collision kernel formulation used below.
All these collision operators are bi-linear in the distribution functions. The first term on the right hand side is due to scattering into the element dv of velocity space and we shall abbreviate it by defining the collision kernel ("redistribution function") C as a proper integral over pre- and post collision velocities of species b-particles:
 δf δt |gain = ⌠⌡ d3v′ C(v′→ v)|v′| f(v′)
(1.3)
The kernel C can be a quite complicated integral, as is involves not only multiple differential cross sections, but also, possibly, particle multiplication factors, e.g. in case of fission by neutron impact, dissociation of molecules by electron impact, or stimulated photon emission from excited atoms. It can also include absorption, in which case the post collision state must be an extra "limbo" state outside the phase-space considered. Due to both particle multiplication and/or absorption the collision kernel C is not normalized to one, generally.
If the distributions fb are assumed to be given, then the kernel C is linear and the expression above becomes a linear integral operator.
The second term on the right hand side is much simpler, because the function f(v) can be taken out of the integral. We even take the product |v| ·f before the integral. The remaining integral is then just the total macroscopic cross section Σt, i.e., the inverse local mean free path (dimension: 1/length). It is solely defined by total cross sections and independent of particle multiplication factors, since we only consider binary collisions (exactly two pre-collision partners always).
This term is then, often, taken on the left hand side of the Boltzmann equation with a positive sign, in the form:
 δf δt |loss = Σt(r,v) |v|f(v)
(1.4)
For the linear case (fb given) this "extinction coefficient" Σt is independent of the dependent variable f = fi0, and this term (out-scattering) just describes the loss of particle flux of i0 particles due to any kind of interaction of them with the host medium.
With these formal substitutions the Boltzmann equation takes a form which is often more convenient, in particular in linear transport theory:

 [ ∂ ∂t + v·∇r + F(r,v,t) m ·∇v] f(r,v,t) +Σt(r,v) |v|f(v)
 =
 ⌠⌡ d3v′ C(v′→ v)|v′|f(v′)
 +
 Q(r,v,t)
(1.1b)
Often the characteristic time constants for neutral transport phenomena are very short (μs), compared to those for plasma transport (ms). We can, therefore, often neglect explicit time dependence in the equations describing the neutral particles. This is done in most applications. The extension to time-dependent problems is rather straight forward and the procedure in the EIRENE code for such cases is described in 1.7 below.
A computationally crucial simplification is the neglect of all neutral - neutral interactions, rendering the equation for f(x) linear. For the status of options in EIRENE to deal with "self-collisions" see below, section 1.9.
For stationary (time-independent) problems the scalar transport flux ("angular flux") Φ, where
 Φ(x) : = | v| ·f( r, v, i) ,
(1.5)
is sometimes used in preference to the distribution f(x) as dependent variable.
In particular for stationary (∂/∂t = 0) and force free (F = 0) problems, as e.g. often encountered in linear transport theory such as neutronics, neutral particle transport in plasmas, etc., the transport equation then reduces to the more compact form:
 ∇r Φ(r,v,t) +Σt(r,v)Φ(r,v) = Q(r,v,t) + ⌠⌡ d3v′ C(r;v′→ v)Φ(r,v′)
(1.1c)
Alternatively, in computational domains with non-vanishing collisionality (i.e., if Σt(x) ≠ 0 everywhere) the (pre-) collision density Ψ is used, i.e.,
 Ψ(x) = Σt (x) ·Φ(x) = νt(x) ·f(x) ,         Σt = Σt(r,E,Ω,i) = νt(r,E,Ω,i)/|v|
(1.6)
where, again, the "macroscopic cross section" Σt is the total inverse local mean free path (dimension: 1/length), and νt is the collision frequency (dimension: 1/time). This cross section can be written as a sum Σt = ∑Σk over macroscopic cross sections for the different types (identified by the index k) of collision processes. Further details about this "macroscopic cross section" and its relation to the conventional ("microscopic") cross sections are given below, section 1.3.5.

1.2  The linear integral equation for the collision density Ψ

By formally integrating the characteristics for (1.1c) the same transport equation can also be written in integral form. This formal procedure is outlined below in paragraph 1.2.1.
The resulting integral equation is often most conveniently written for the collision density Ψ (1.2) rather than for transport flux Φ:
 Ψ(x) = S(x) + ⌠⌡ dx ′Ψ(x ′) ·K(x ′→ x) .
(1.1d)
This equation has the general form of the backward integral equation of a Markovian jump-process and it is therefore particularly well suited for a Monte Carlo method of solution. The formal relation between the integro-differential form (1.1c) and this integral form is very useful to generalize the Monte Carlo procedure, e.g., to time-dependent equations, and to Boltzmann-Fokker-Planck equations (which contain diffusive contributions or diffusive approximations for some processes, in addition to the jump processes described by the Boltzmann collision integral). It also allows to make connection to the so called "Greens-functions Monte Carlo" concept (originally developed for quantum mechanical problems involving solutions to the Schrödinger equation). A corresponding discussion is postponed to section 1.2.1. A direct intuitive interpretation of the integral equation is already sufficient to understand the Monte Carlo method of solution and shall be given first.
In (1.1d) x′ and x are the states at two successive collisions (jumps). The integral ∫dx′ in (1.1d) is to be understood as an integral over all initial coordinates, i.e. over the entire physical space, the full velocity space and a summation over all species indices. The transition kernel K is usually decomposed, in our context, into a collision- and a transport kernel, i.e., C and T , where
 K(r′, v′, i ′→r, v, i) = C( r′;v′, i ′→ v, i) ·T(v, i; r′→ r) .
(2.7)
The kernel C is (excluding normalization) the conditional distribution for new co-ordinates (v,i) given that a particle of species i′ and with velocity v′ has undergone a collision at position r′. This kernel can further be decomposed into:
 C( r′, v′, i ′→v, i) = ∑ k pk  Ck ( r′; v′, i ′→ v, i)  ,  pk = Σk Σt
(2.8)
with summation over the index k for the different types of collision processes under consideration and pk defined as the (conditional) probability for a collision to be of type k. The normalizing factor
 ck (x ′) = ∑ i ⌠⌡ d3v  Ck (r′, v′, i ′→v, i )  ,  Ck = 1 ck Ck
(2.9)
gives the mean number of secondaries for this collision process. The function Ck then is a conditional probability density. The particle absorption process can conveniently be described by adding an "absorbing state" xa to the μ-space (generally referred to as "one-point compactification of this space" in the language of mathematical topology). This "limbo state", once it is reached, is never left again if the kernels T or C are employed as transition probabilities. Formally, an additional collision kernel Ca (xxa ) and an absorption probability pa = Σa / Σt must be included in the collision kernel. The quantity Σa comprises all collision processes with no next generation particles.
The kernel T describes the motion of the test particles between the collision events. Let, again, Ω′ denote the unit vector in the direction of particle flight v′/|v′|, and let Ω2′ and Ω3′ be two further unit vectors such that these three vectors form an orthonormal basis at the point r′. Neither velocity nor species change along the transition described by T, i.e. v′=v and i′−i. Omitting the corresponding delta functions in velocity space and the Kronecker delta δii the transport kernel T then reads as follows:
 (T(l)=)  T ( v′, i′; r′→ r )
 =
 Σt ( r,v, i) ·exp ⎡⎣ − ⌠⌡ l=Ω′( r−r′) 0 ds Σt ( r′+ sΩ′, v′, i′) ⎤⎦
 ·δ(Ω2′( r−r′)) ·δ( Ω3′(r− r′)) ·H (Ω′( r− r′))
(2.10)
 =
 Σt ( r,v, i) ·F( v′, i′; r′→ r )     0 ≤ l ≤ ∞
(2.11)
with H(x) = 0 if x ≤ 0 , and H(x) = 1 if x > 0, the Heaviside step function. The two remaining delta functions restrict the motion to a path in the direction of the initial velocity v′.
Thus, although strictly being a conditional (on x′) distribution in phase space, for an infinite medium T can be interpreted as the distribution density for the distance l for a free flight starting from r′ to the next point of collision r = r′+ l ·Ω′.
For a finite medium this distribution can be generalized to (writing shorter l for (r′+lΩ′,v′,i′):
 T(l)
 =

 Σt ( l)
 ·exp ⎡⎣ − ⌠⌡ l 0 ds Σt ( s ) ⎤⎦ ,
 l < lmax
 δ(l−lmax)
 ·exp ⎡⎣ − ⌠⌡ lmax 0 ds Σt ( s ) ⎤⎦ ,
 l ≥ lmax
(2.12)
Here lmax denotes the distance along the flight direction from r′ to the boundary for the computational domain, to any internal surface at which the test flight shall be stopped, e.g. for scoring surface fluxes there, or even to the cell boundary. The "Transport Kernel" T has dimension: [1/Dimension of phase space] (T(x′→ x)dx is a probability). The function F (Dimension: [length times Dimension of T]) defined in expression (2.9b), and analogously from (2.10) for a finite medium, will turn out to be the relevant Green function for the transport problem, see Section 1.2.1. The integral
 α(r′,r) = ⌠⌡ Ω( r− r′) 0 ds Σt ( r′+ s Ω)
in equation (2.9a) is well known as "optical thickness of the medium" in linear transport theory.
The inhomogeneity S in equation (1.1d) is, excluding normalization, the distribution density of first collisions, whereas the integral term in equation (1.1d) describes the contribution to Ψ from all higher generations. The quantity S can be written as:
 S(x) = ⌠⌡ dx ′Q(x ′) ·T(x ′→ x) ,
(1.6a)
with a source density Q. As the problem is linear, Q can be normalized to 1 and, thus, Q can be considered a distribution density in μ -space for the "primary" birth points of particles, as, e.g., opposed to the "secondary" birth point distribution (or "post collision density") χ of particles after a collision event

 χ(x) = ⌠⌡ dx ′Ψ(x ′) ·C(x ′→ x).
(1.6b)
It can be shown that a unique solution Ψ(x) exists subject to appropriate boundary conditions and under only mild restrictions (basically on the constants ck and pa ) to ensure that the particle generation process stays sub-critical.
Usually, a detailed knowledge of Φ or Ψ is not required, but only a set of "responses", R, defined by
 R =     < Ψ| gc >     = ⌠⌡ dx Ψ(x) ·gc (x) ⎛⎝ =     < Φ| gt >     = ⌠⌡ dx Φ(x) ·gt (x) ⎞⎠ ,
(2.13)
where gc (x), gt (x) are given "detector functions".
For example all terms in the plasma fluid equations resulting from neutral plasma interaction can be written in this way. This can be seen by considering a numerical grid, composed of M mesh cells (spatial and/or temporal), for the numerical solution of the fluid equations. The detector functions for many responses needed for fusion plasma applications are hard wired in EIRENE, generalization to any arbitrary response by resorting to "user defined detector functions" is described in Section 3.2.1 for tracklength estimates, and in Section 3.2.2 for collision estimates.
Lets therefore define an entire set of detector functions gm, one for each mesh cell of an external code, each including a characteristic function
 gm=g ×chm(r,t) ,     m=1,2,...,M,
(2.14)
i.e., chm (r,t)=1 inside the numerical mesh cell (or time interval) labeled with the cell index m, and chm (r,t) = 0 outside this cell. Thus profiles of cell volume averaged responses are readily obtained in a single Monte Carlo run.
Estimates of surface fluxes, or point estimates e.g. in time, are also included in this concept if proper use is made of delta functions to reduce dimensionality of the response:
 gm,α=g ×chm(r,t) ×δα(r,t) ,     m=1,2,...,M,
(2.15)
Here, e.g., a surface cell m would discretise a surface S characterized by a surface delta function δS.

Monte Carlo estimates of the volumetric source terms in the fluid equations due to the trace particles (neutral particles, but also trace ions described by the Monte Carlo procedure) are therefore cell averages, surface averages, time averages or point estimates (e.g. in time), averaged over a sub-manifold with reduced dimensionality, however, only if the probability for particle histories crossing this sub-manifold is not zero. Surface averages in EIRENE are described in Section 3.2.4, point averages in time are described in Section 3.2.3. This option to reduce dimensionality of responses usually excludes point estimates in real space, for which special estimators not mentioned here would be required.
Nevertheless, in concluding, depending upon the numerical algorithm in the fluid code, one may then have to interpret these Monte Carlo estimates properly, e.g. they may have to be interpolated to the grid points, or be properly rescaled in case of time-dependent applications.

1.2.1  The Green function concept

We return to Equation (2.9b), where the function F(v,i;r′→ r)) was defined. For a finite domain, see (2.10), and introducing again the shortcut l for (r′+lΩ,v,i) this becomes:
F(l) =

 exp ⎡⎣ − ⌠⌡ l 0 ds Σt ( s ) ⎤⎦ ,
 l < lmax
 0
 l ≥ lmax
(2.16)
to be written

1.3  Monte Carlo solution of equation 1

A statistical solution to equation (1.1b) is straight forward, because it is formulated in probabilistic terms as follows.
A discrete Markoff chain is defined using Q as an initial distribution and L = T ·C (order of C and T reversed compared to K) as a transition probability. Histories ωn from this stochastic process are generated according to ωn = (x0, x1, x2, ... , xn), (where xj = xa for all jn and xixa for all i < n), with xn being the first state after transition into the absorbing state xa. x0 denotes the initial state distributed as described by Q. Thus, the length n of the chain ωn is itself a random variable. A random sampling procedure to generate such chains is carried out in Monte Carlo codes by converting machine generated (pseudo-) random numbers ξi1 , ξi2, ... into random numbers with the distributions Q, T and C. Having computed N (several thousand) chains ωi, i=1,2,...,N, the responses R are estimated as the arithmetic mean of functions ("statistics", or "estimators") X(ω), i.e.
 R ≈ ~ R = 1 N N∑ i=1 X(ωi) .
(3.17)

1.3.1  Unbiased estimators

One possible choice for X(ω) is the so called "collision estimator" Xc,
 Xcg ( ωni) = n∑ l=1 gc (xl ) · l−1∏ j=1 c(xj ) (1 − pa(xj)) .
(3.18)
This estimator is, for example, used in the DEGAS code [10]. It can be shown that the statistical expectation E(Xc) produces:
 R = E(Xc ) = ⌠⌡ d( ω) Xcg ( ω) h( ω)
(3.19)
with h(ω) being the probability density for finding a chain ω from the Markoff process defined above, i.e. Xc is an unbiased estimator for R.
Other estimators ("track-length type estimators") are employed frequently. These estimators are unbiased as well but have higher moments different from those of Xc. Instead of evaluating the detector function gc(x) at the points of collisions, xl , (as Xc does,) they involve line integrals of gt(x) along the trajectories, e.g.,
 Xtg ( ωni ) = n−1∑ l=0 ⎧⎨ ⎩ ⌠⌡ xl+1 xl ds  gt (s) ⎫⎬ ⎭ · l−1∏ j=1 c(xj ) (1 − pa(xj)) ,
(3.20)
with R = E(Xt ) = E(Xc). It has been shown, e.g., in [7], that the collision estimator, derived not for the pre- but for the post-collision density integral equation results in a track-length type "conditional expectation estimator" Xe, which, together with Xc and the "traditional" track-length estimator Xt, may be used as one further option in the EIRENE code.
This estimator Xe is obtained from Xt by extending the line integration, which is restricted to the path from xl to xl+1 in formula (3.17), to the line segment from xl to xend. Here xend is the nearest point on a boundary along the test flight originating in xl. I.e., the line integration may be extended into a region beyond the next point of collision, into which the generated history would not necessarily reach. This "conditional expectation estimator" reads:
 Xeg ( ωni ) = n−1∑ l=0 ⎧⎨ ⎩ ⌠⌡ xend xl ds  gt (s) ·exp ⎛⎝ − ⌠⌡ s 0 ds ′Σt ( s ′) ⎞⎠ ⎫⎬ ⎭ · l−1∏ j=1 c(xj ) (1 − pa(xj)) ,
(3.21)
If xend is taken to be the nearest point on each mesh cell boundary, then the estimator (3.18) reduces, as a special case, to the method employed by the NIMBUS code [11]. However, the estimator (3.18) is more general, as the integration may be extended over arbitrarily many cells. The length of the integration path is controlled in EIRENE by an input flag WMINC, see section 2.10.

1.3.2  Scaling of tallies

With increasing N, the number of Monte Carlo histories, the unbiased estimators X given in previous section 1.3.1 provide arbitrarily precise approximations for the responses Rg=〈Ψ|g〉, for detector function g and dependent variable Ψ (equation 1.1d). The precise meaning and physical dimensions of g and Ψ depend upon the problem at hand, e.g. also on symmetry (ignorable spatial coordinates, or time).
Because the inhomogeneous part of the transport equation (1.1d) S was assumed to be normalized to one (for source sampling),
 S (r,v,t)= 1 s S(r,v,t),     s = ⌠⌡ d3rd3vdt S(r,v,t)
(3.22)
this normalization factor s has to be multiplied to estimators (3.14) to turn responses Rg (2.11) into estimates with correct dimensionality and in absolute units.
Mostly the responses of interest are intensive quantities, such as density, pressure, collision density (number of collision per unit time and volume) rather than the extensive quantities obtained by the phase space integration in (2.11), such as volume, total energy or number of particles. Thus by dividing the estimate by the appropriate cell-volume Vm (in space-time), see (2.12), finally the profiles of cell averaged intensive quantities (e.g. density profiles) are obtained. The final unbiased estimate, in absolute units, for any of the unbiased estimators Xg for detector function g discussed above, then is:
 Rg = 〈Ψ|g〉 ≈ ~ R g (N) = s Vm×N N∑ i=1 Xg(ωi) ,     N >> 1
(3.23)
with N being the number of Monte Carlo histories. Correct interpretation of the results of an EIRENE run hence requires knowledge of the ratio s/Vm. The total source strength scaling factor s, i.e. the integral of the inhomogeneous part S of the governing Fredholm integral equation, is input to EIRENE, variable "FLUX" in input block 7 (see 2.7). The space-time cell volumes Vm are either automatically calculated by EIRENE in subroutine "VOLUME", for cells belonging to standard grid (input block 2, Section 2.2), or are provided from external considerations by a number of options, e.g. input blocks 3b, 5, or 8. In this latter case great care is needed, however, that the cell volumes are proportional to the volumes as seen by the test flights for otherwise not only unscaled profiles, but even wrong (biased) profile shapes will be obtained.

Scaling in problems with ignorable coordinates

As seen in the previous paragraph, the absolute values of estimates are scaled with the ratio s/V of source strength to cell volume. Depending on ignorable coordinates in any particular problem, or on whether stationary (i.e.: time is an ignorable coordinate) or time dependent transport problems are considered , the dependent variable Ψ in the governing integral equation (1.1d), the source strength s and "cell volume" V may have different interpretations:
stationary (time independent) problems
one ignorable spatial coordinate   In stationary problems with one ignorable coordinate, say, the z-coordinate, the source strength s (input flag "FLUX") is the flux, particles per unit length dz in direction of z. Likewise, the volume V is per unit (mostly "toroidal") length, i.e. if dz=1 then V is the cell area in the two remaining coordinates x,y. In

1.3.3  Efficiency (FOM)

The efficiency for Monte Carlo Codes is the inverse of the figure of merit (FOM) of the calculation, defined as
 FOM = statistical  variance ·computing  cost.
(3.24)
Note that FOM should be approximately independent of the running time, because the number of histories generated is (excluding overhead) proportional to the CPU costs and inversely proportional to the statistical variance σ2.
It is one of the major advantages of Monte Carlo methods over other numerical schemes that the error estimates, empirical variances [(σ)\tilde]2, are directly provided by the method itself, not requiring any further considerations. The options to activate evaluation of statistical variance in an EIRENE run for any computed quantity (tally) are described in input block 9 (section 2.9). EIRENE provides numerical and graphical output for the "empirical relative standard deviation", in %.
 ~ σ g,rel (N) = ~ σ g (N)/ ~ R g (N) ×100
(3.25)
and [(R)\tilde]gN is the Monte Carlo estimate for tally Rg based on N Monte Carlo histories (3.20). Here Rg may be any volume or surface averaged tally for detector function g, either tracklength- collision- or snapshot estimated. The variance σ21 per history is obtained as the (unbissed) estimate ("empirical variance") as:
 σ21 ≈ ~ σ 21 (N) = 1 N−1 N∑ i=1 (Xi− X )2= 1 N−1 ⎡⎣ N∑ i=1 Xi2 − 1 N ⎛⎝ N∑ i=1 Xi ⎞⎠ 2 ⎤⎦
(3.26)
where Xi = Xgi) is the contribution of Monte Carlo history ωi to estimator Xg for tally g, and [ˉ(X)] is the arithmetic mean [ˉ(X)]=1/NXi over all histories. N is the number of (statistically independent) Monte Carlo histories. The subscript g is omitted here and from now on. The variance per history is turned into the final estimate for the variance by the "law of large numbers" (and the "central limit theorem" of probablity theory):
 ~ σ 2 (N)= 1 N ~ σ 21 (N)
(3.27)
For large N the variance per history [(σ)\tilde]12 converges to a constant (namely to σ21), and the final Monte Carlo estimate of variance [(σ)\tilde]2 (3.24) has the expected probabilistic 1/N scaling, i.e. the empirical standard deviation scales as: [(σ)\tilde]  ∼ 1/√N.
Multiple cell crossings   One should note that except in simple 1D cases in general one history ωi can contribute more than once to the estimate for a particular cell of the grid. E.g. test particles can cross one cell m (see (2.12b) more than once, with each cell crossing "j" leaving a score Xi,j. Then: Xi=∑j Xi,j. The summation in (3.14) and hence also the second sum on right hand side of (3.23) can still be accumulated "on the fly" while generating the Monte Carlo histories.

Sampling, Non-analog sampling

One often encounters in literature the description of special, "new" Monte Carlo techniques, that are greatly superior to so called "standard methods".
It is true that one can devise very intelligent methods to optimize performance, but each optimization always only works well for one very particular problem, it can equally well entirely wreck the performance of an only slightly different case.
For a general purpose Monte Carlo solver for transfer problems, as EIRENE, therefore, no general advise can be given, although the performance often could greatly be improved by adapting the method to a particular problem. In particular, non of these "intelligent methods" have been (and will ever be) hard wired into EIRENE.
EIRENE contains a set of such methods, referred to as "non-analog" methods and controlled by the flags in input block 9, see section 2.9. Activating those must be accompanied by a very careful statistical analysis of results, not just the run time per particle, nor even the standard deviation alone, suffices.
We first note here that any random distribution function f(x) arising during the course of generating the chains ω ("particle histories") can be replaced by another "non-analog" function, g(x), if a weighting factor,
 w(y) = f(y) / g(y)
(3.28)
is included in the formula for the estimator X(ω), y being an actual random number generated from g . This choice can in some cases increase the efficiency of the algorithm.
The only restrictions on the choice of non-analog functions (besides practical ones) are:

 if     g(x) = 0,        then     f(x) = 0;
(3.29)
and conditions to ensure that the non-analog process remains sub-critical as well.
The condition 3.26 is checked in an EIRENE run whenever non-analog distributions are applied, and, if violated, an error exit with the message: "violation of Radon-Nikodym condition" occurs.
Note: a violation of condition 3.26 can not be detected otherwise, e.g., by monitoring overflow in the weighting factor w. Of course, values y resulting in a zero in the denominator of w (equation 3.25) have sampling probability zero according to distribution g(x), and, hence, are never sampled. Still the results would be biased.
Note: when inspecting an EIRENE geometry plot with particle trajectories (or even a movie: NLMOVIE=TRUE, input block 1) in order to get an intuitive feeling for the particular transfer process considered, then all non-analog sampling options must first be turned off, (NLANA flag in input block 1), for otherwise the pictures may be grossly misleading.

We consider the procedures for random sampling from univariate distributions as known, and refer to the many textbooks on that, in particular to the "random sampling library" [12]. If none of the direct methods apply, then still either the non-analog method mentioned above, or the "rejection method" can be used.
Sampling from a multivariate distribution f(x1,x2,...,xn) can always be reduced to a sequence of samplings from univariate distributions, by noting that:
 f(x1,x2,...,xn)=f1(x1) ·f2(x2|x1) ·f3(x3|x1,x2) ·...
(3.30)
Here, f1 is the marginal distribution obtained from f by integrating over all but the first independent variables. It is a univariate distribution of x1.
f2 is a conditional marginal distribution obtained from f firstly by integrating f over all but the first 2 variables x1, x2, and secondly then taking the conditional distribution, conditional on x1 (i.e., the univariate distribution of x2 for given values of x1).
Likewise, f3 is a conditional marginal distribution, n-3 fold integration of f, and then distribution conditional on x1, x2, and so on.
Random sampling from f(x1,x2,...,x3) then proceeds by sampling first x1 from the (univariate) first factor in equation 3.27, then x2 by sampling from the (univariate) second factor, and so on.
One particular example of this scheme is the "TRIM database surface reflection model" (see below, section 1.4.2). There, essentially, tables of the conditional marginal distributions mentioned here are pre-computed with the TRIM code. For convenient random sampling, these tables are stored for the inverted cumulative distribution, i.e., as conditional quantile functions.

Stratified Source Sampling

The EIRENE code resorts to a "stratified sampling" technique. This means that the primary source distribution (see Equation 1.1a) can be decomposed into a sum of independent sources:
 Q = ∑ i Qi
and the solution is obtained by linear superposition of the solutions for each sub-source ("stratum") Qi.
The stratified sampling technique is a well known statistical procedure, see any textbook on statistics or Monte-Carlo methods.
It sometimes can significantly affect the efficiency of a Monte-Carlo run. This can go both ways, depending upon the particular problem and the particular stratification used. Therefore, as with the non-analog options, no general recommendation can be made.
The effects of stratification, however, are usually more easy to assess (predict) than with general non-analog methods.
In order to make connection with the textbooks we need to note here only the following:
EIRENE provides estimates of linear functionals ("moments", "responses")

 < g,φQ > = < Q,φ*g >
(3.31)
As described above, g is the detector function ("weighting function") φQ is the solution to the kinetic transfer equation (Boltzmann eq.) in which Q was used as (inhomogeneous) source term. φ*g is the solution of the corresponding adjoint equation, with g as source term.
The above equation means that EIRENE estimates can simply be regarded as Monte-Carlo estimates of multi-dimensional integrals, namely of the adjoint function φ*g integrated over phase space with Q defining a distribution (measure) in phase space "μQ".
It is then obvious that one can split that integral into a sum of integrals by arbitrarily decomposing the domain of integration, i.e., by decomposing the source Q into sub-sources ("strata") Qi with
 Rg = ⌠⌡ dxφQ ·g = < Q,φ*g > = ∑ i < Qi,φ*g > = ∑ i ⌠⌡ dμQi φ*g
(3.32)
This is the concept of stratified sampling, which, for Monte Carlo particle transfer procedures in particular actually turns out to be a "stratified source sampling".
Further details about this in our particular context are given in the description of input block 7, section 2.7 below.

1.3.4  Source sampling

to be written:
and see special section on "plasma surface recycling source sampling".

1.3.5  Sampling from the transport kernel T

The conditional sampling distribution T(r′→ r) for the next point of collision (event), given a test flight is at position r, traveling with velocity v, is given by (2.9a). Let us introduce the distance of the flight l in direction Ω =v/|v|, hence r =r′+lΩ, then the sampling distribution for l is simply
 T(l)=Σt(r′+lΩ)exp ⎡⎣ − ⌠⌡ l 0 ds Σt(r′+sΩ) ⎤⎦ ,     Σt=1/λt
(3.33)
Distances l are hence readily sampled from this exponential distribution using the inversion method and the cumulative distribution GT of T
 ξ = GT(l)=1−F(l)=1−exp ⎡⎣ − ⌠⌡ l 0 ds Σt(r′+sΩ) ⎤⎦
(3.34)
with ξ (and hence also 1−ξ) being a uniformly distributed random variable. Therefore, cumulating exp[−∫0l ds Σt(r′+sΩ)] along a flight and comparing it with ln(ξ) with a uniform random number ξ provides a random sample of the free flight distance l.
Sampling the distance l of a free flight, for test particle species "i" with velocity v, from the transport kernel T requires knowledge of the ("monochromatic") mean free path λ = λ(r,|v|) along the trajectory, and, hence, of the "macroscopic cross section" Σ = Σ(r,|v|), which also may vary along the trajectory.
The "macroscopic cross sections" Σ (by abuse of language, dimension: 1/length) with or without subscript k (the label for the type of the collision process) or subscript t (for "total", the sum over k), are defined as:
 Σ = 1 λ = ν |v| ,       λ = |v| ν
(3.35)
with the mean free path λ and the collision frequency ν:
 ν(r,v,i) = nb 〈σ·vrel〉b,        vrel = |v−vb|
(3.36)
Here nb and vb are the target particle density and velocity, respectively, and the brackets 〈...〉b denote averaging with the velocity distribution fb(r,vb) of the background medium species b.
σ is the total scattering cross section, σ = σ(|vrel|).
Hence, e.g., for all isotropic target distributions fb one has:
 ν = ν(r,E,\not Ω,i)       and        Σ = Σ(r,E,\notΩ,i).
Furthermore, let fb be a drifting Maxwellian, with temperature Tb and drift Vb: fb = fb(Tb,Vb).
Then:
ν = ν(Tb,|vVb|,i)
Finally, if the thermal speed vbtherm in this Maxwellian: vbtherm = √{Tb/mb} is very large compared to typical values of |vVb| (i.e.: nearly incompressible flow conditions), then:
ν ≈ ν(Tb,0,i) = nb 〈σ·vb
and here 〈σ·vb is the usual "Maxwellian rate coefficient" for the collision process, taken at temperature Tb of the background particles (species b) and at zero velocity of the test particles (species i). For example, in case of neutral atom or molecule collisions with electrons this approximation is usually valid.

1.3.6  Sampling from the collision kernel C

The general rules for sampling from multivariate distributions outlined in section 1.3.3, equation 3.27 also apply here.
It is not necessary, but usually most convenient and intuitively most clear, to decompose the kernel C according to the physical nature of the various collision processes taken into account.
The collision kernel C was written as (see equation 2.7)
 C( r′, v′, i ′→v, i) = ∑ k pk Ck ( r′; v′, i ′→ v, i)  ,  pk = Σk Σt
with summation over the index k for the different types of collision processes under consideration and pk defined as the probability for a collision to be of type k. The normalizing factor
 ck (x ′) = ∑ i ⌠⌡ d v Ck (r′, v′, i ′→v, i )  ,  Ck = 1 ck Ck
was the mean number of secondaries for this collision process. The function Ck then is a conditional probability density for (v,i), given the pre-collision coordinates (r′,v′,i′).
Sampling from C (given that a collision of some kind has occurred) proceeds most conveniently by firstly sampling the type k of the process from the discrete distribution
 (pk,     k=1,..,K),
and then by sampling the post-collision state of the test-particle from Ck. Finally the weight of the test-flight after collision is increased (or reduced) by multiplying it with the normalization constant ck.
The decomposition of C into sub-kernels Ck is somewhat arbitrary, and the criteria may also be computational aspects rather than the different physical nature of the various processes taken into consideration.
Sampling (v,i) from Ck can be done along the same idea: first factoring out the discrete distribution for the species index i, and then by sampling v from Ck,i, the conditional distribution for v given the type of the process is k and the species of the post collision trajectory is i. Let (nk,i,i′=1,...,Ik) be the number of post collision particles of species i′ from process k, and nk = ∑i nk,i. Then
 ( ~ n k,i′ i′=1,...Ik) = (nk,i′/nk,     i′=1,...Ik)
is the discrete distribution for the post-collision particle species i′, given a collision process of type k.

1.4  Surface Reflection Models

Within the terminology fixed in section 1.2 the effects of the interaction of neutral particles (and escaping ions) with surfaces surrounding or inside the computational volume can be described by additional boundary conditions. It is, however, much more convenient for Monte Carlo applications to describe this interaction as just one special type of collision process in the collision kernel C of equation (2.7), i.e. as
 C = pw ·Cw + ∑ k pk Ck ,
(4.37)
with pw (r) = 1, if r is a point on a surface, and pw (r) = 0 elsewhere. The transport kernel T can be modified such that the maximum length lmax for free flight to the next point of collision is the distance to the nearest surface along the track of the particle. The reflection kernel Cw is further decomposed it into a kernel Cwf (fast particle reflection), Cwt (thermal particle re-emission) and Cwa (particle absorption) with respective probabilities, pf , pt and pa, such that
 Cw = Cwf + Cwt + Cwa = pf Cwf + pt Cwt + paCwa ,     pf + pt + pa = 1 .
(4.38)
Here C are, as above, the normalized versions of the collision kernels C.
We discuss in detail only the fast particle reflection model Cwf since Cwt is generally a mono-energetic and cosine angular distribution or a Maxwellian flux distribution at wall temperature, and Cwa only describes the transition into the "limbo state" xa "absorbed particle".
It is usual in neutral gas models to replace Cwf with numerical or analytical fits to moments of Cwf, e.g., such as the particle and energy reflection coefficients CP and CE. The particle - wall interaction is often supplemented by simple assumptions on the angular distribution (cosine or specular) of the particles leaving the surface.
One has for the particle reflection coefficient:
 CP = ∑ i ⌠⌡ d v Cwf (r; v′, i ′→ v ,i)
(4.39)
Note that CP = pf = cwf in the terminology of equation 2.8.
The energy reflection coefficient is defined as:
 CE = 1 E ′ · ∑ i ⌠⌡ d v E ·Cwf (r; v′, i ′→ v,i) = 1 E ′ · ~ C E
(4.40)
Here E = m v2 / 2 is the energy of the reflected particle. Hence CE = [ˉ(E)]/E′·CP, with [ˉ(E)] denoting the mean energy of the reflected particles.

1.4.1  The Behrisch matrix reflection model

One model, the matrix model due to Behrisch [13], is hardwired in EIRENE (and in many other codes such as AURORA, BALDUR, TRANSP). This model has been used as a standard option in many benchmarks [14]. The model is based upon a table of values for CP (E ′) and a stochastic matrix B for the transition E ′→ E, for H atoms striking a stainless steel target at normal incidence. For other incident angles ϑin the moments are modified in EIRENE by the relations:
 CP ( ϑin ) = 1 − ( 1 − CP⊥) ·cose1 ( ϑin )
(4.41)
and
 CE ( ϑin ) = 1 − ( 1 − CE⊥) ·cose2 ( ϑin ) .
(4.42)
This latter moment equation can be satisfied by modifying the tabulated stochastic matrix B( E ′→ E) in many different ways. In EIRENE, we first sample E from B (given E ′) and then set
 E = E ′− ( E ′− E ) ·cose2 ( ϑin ) .
(4.43)
The linearity of E in E gives an expectation value [ˉ(E)] with [ˉ(E)] / E ′·CP( ϑin ) = CEin ), with CP and CE as given by equations (4.38) and (4.39), respectively.
To describe the angular distribution of re-emitted particles, let ex , ey , ez define an orthonormal basis at the strike-point of a test-particle on a surface. Here let ex be the unit vector parallel to the outer surface normal and let the incident particle travel in the local xz-plane. Furthermore, let ϑr and ϕr be the polar and azimuthal angles sampled from a cosine distribution around − ex. The speed unit vector Ω of the reflected particle in EIRENE is then given by
 Ωx
 =
 − cos( ^ ϑ r ) f1 ( ϑin) + sin( ^ ϑ r )sin( ϕr) f2 ( ϑin ) ,
 Ωy
 =
 sin( ^ ϑ r )cos( ϕr ) ,
 Ωz
 =
 sin( ^ ϑ r ) sin( ϕr )f1( ϑin ) + cos( ^ ϑ r ) f2(ϑin ) ,
(4.44)
with
 sin( ^ ϑ r )
 =
 f1( ϑin) sin( ϑr ) ,
 f1( ϑin )
 =
 cose3 ( ϑin ) ,
 f2( ϑin )
 =
 ( 1 − f1 ·f1 )[1/2] .
In total we have defined three parameters e1, e2, e3. The first two parameters e1 and e2 account for particle and energy reflection coefficients that continuously increase with incident angle (recommended: e1=1, e2=0.5 from comparison with TRIM code calculations, see below).
The third parameter e3 can be used to take into account an increasing contribution of specular reflection at near grazing incidence. e3=0 gives a pure cosine distribution, e3=1 reproduces the distribution given in [10]. In general this model describes a fraction of specular versus cosine reflection which increases with ϑin, and in the (purely mathematical) limit of ϑin = 90°, pure specular reflection. Furthermore, the parameter e3 controls the speed at which the specular part aspec (Ω) of the angular distribution a(Ω) = aspec ( Ω) + acosin ( Ω) increases with ϑin .
The interaction of neutrals and ions with a solid surface for target-projectile combinations other than H onto stainless steel are modelled using a reduced energy scaling.

1.4.2  TRIM code database reflection models

Complementary to the surface model described above, EIRENE is coupled to computer generated reflection databases, produced with the TRIM code (see references [15] and [16]). At present a discretized form of the reflection kernel Cwf, including all correlations between the vector components of v′ and v, is available for the target-projectile combinations listed in table 1.1. ZNML is a target material flag set for each surface in the input file, blocks 3a and 3b, and is described below in block 6b, "Data for local reflection models". These databases are, essentially, tables of conditional quantile functions obtained from large random samples of test particle trajectories in the solid.
Table 1.1: Surface reflection database (TRIM Code)
 No Projectile Target ZNML 1 H Fe 5626. 2 D Fe 5626. 3 H C 1206. 4 D C 1206. 5 He Fe 5626. 6 He C 1206. 7 T Fe 5626. 8 T C 1206. 9 D W 18474. 10 He W 18474. 11 H W 18474. 12 T W 18474. 13 D Be 904. 14 D Mo 9642. 15 Ne C 1206. 16 Ne Be 904. 17 H Cu 6429. 18 H Mo 9642. 19 T Mo 9642. 20 He Mo 9642.
We consider these "database -" reflection models to be the most complete option available at present and recommend to further extend the databases to additional projectile - target combinations, and in particular to sputtering and self-sputtering data (such as, e.g., C onto C, or Be, etc.).
For an INTOR benchmark case for neutral particle Monte Carlo Codes [14] the choice of a particular surface reflection model had no strong impact on the results. This was shown by running the EIRENE code with all three models ("Behrisch Matrix", TRIM database and MARLOWE database (at that time available as a third option. Meanwhile the MARLOWE database is not available anymore).
It was found that this had little influence on global parameters as, e.g., the pumping efficiency for the particular geometry (5.4 %, 5.3 % and 5.2 % respectively.). On the other hand, it can easily be imagined that in other cases (different plasma conditions, more detailed geometries) the choice of the reflection kernel Cwf can be more important for the result. Certainly in all cases, in which the results are sensitive to the neutral particle velocity distribution near surfaces, details of the reflection kernel matter. This is, for example, the case for computed line shapes of certain electronic transitions (Hα lines etc.), see e.g., reference [17].

1.5  Surface sources

It has been found as a result of benchmark calculations with several different neutral gas transport Monte Carlo codes on the same input data that the main divergence in the results from the different codes tested was due to different primary source terms Q in the transport equation [equations (1.1a,1.1b) and (1.6a)] , although target fluxes and temperatures had been prescribed to be identical. Here we describe one particular EIRENE surface source option in detail. The way the source is modelled by EIRENE allows consistency with kinetic boundary conditions for plasma fluid equations at plasma facing surfaces. To describe this, we first write the volumetric source Q (particles per unit time and volume) in terms of a surface (wall) flux Γw and surface delta-function:
 Q( r, v, i ) = Γw (r, v, i ) ·δ( π) ,
(5.45)
where π is a coordinate normal to the surface Sw: π(r)=0 , i.e.: eπ = ∇rπ(r), pointing away from the computational volume. Γw is a surface source (flux-) density distribution, backward directed (into the computational volume, i.e., only vπ < 0).
A typical example is a sonic parallel plasma flow (Mach number M||=1 along the magnetic field lines) entering the magnetic presheath. In this magnetic presheath the ion trajectories are bend over such that at the entrance of the electrostatic sheath the sound speed is already reached with respect to the surface normal: Mπ=1.
To convert from a forward directed ion flux distribution ΓSh,ion at the plasma-sheath interface SSh to a neutral flux distribution back from the wall surface Sw, we fold ΓSh,ion with a (magnetic plus electrostatic) sheath transmission kernel TSh and with the surface reflection kernel Cw, equation 4.34:
 Γw ( r, v, i )
 =
 ∑ i ′ ⌠⌡ ⌠⌡ ⌠⌡ d3 r′d3 v′d3v" ΓSh,ion ( r′,v" , i ′)
 ·TSh (i ′; r′,v" → r, v′) ·Cw ( r;v′, i ′→ v, i )
(5.46)
The summation over i ′ is performed over all ion species entering the sheath region.
A typical example for ΓSh,ion is a sonic parallel plasma flow (Mach number M||=1 along the magnetic field lines) entering the magnetic pre-sheath. In this magnetic pre-sheath the ion trajectories are bend over such that at the entrance of the electrostatic sheath the sound speed is already reached with respect to the surface normal: Mπ=1.
A commonly adopted simplification is made by ignoring the typically small spatial volume between the plasma sheath interface and the wall surface. Hence: TSh factorizes into a 3D delta function δ3(rr′) and a deterministic transition kernel for the velocities v" → v′.
Random numbers from Γw are generated conveniently by sampling, firstly, (r′, v" , i ′) from ΓSh,ion, next generating (r,v′) from TSh and then, finally, (v, i) from the conditional distribution Cw .
A frequently employed assumption for ΓSh,ion is a forward drifting Maxwellian ion flux distribution at the interface between plasma and electrostatic sheath SSh. In this case the bending of trajectories from parallel to normal (to the surface) sonic flow in the magnetic pre-sheath is already included in ΓSh,ion , e.g. via boundary conditions in a CFD model for plasma ion flow, and the action of operator TSh is solely the acceleration in the electrostatic sheath. Such a forward drifting Maxwelllian ion flux distribution is given as
 Γ0Sh,ion = 1 α2 ·vπ · exp ⎛⎝ − 1 2 β2 (v− Vd )2 ⎞⎠ ,    vπ > 0
(5.47)
with β = √{ [(kT)/(m)] }, normalizing constant α2 and with the drift velocity Vd. The superscript 0 indicates normalization of the flux to one.
Consider the three orthogonal velocity components v1,v2,vπ, then Γ0Sh,ion(v) factorizes in the three uni-variate distributions
 Γ0Sh,ion(v)=Γ0Sh,ion,1(v1)·Γ0Sh,ion,2(v2) ·Γ0Sh,ion,π(vπ)
(5.48)
The first two are simply shifted Gaussians. However, sampling the π-component vπ of v from ΓSh,ion is not trivial because the cumulative distribution function F of a 1D forward drifting Maxwellian (with cut-off at vπ=0):
 F ( vπ ) = ⌠⌡ vπ 0 dvπ ′ Γ0Sh,ion ( vπ ′)
(5.49)
cannot be inverted explicitly. Instead we sample this component [(v)\tilde] π from a "non-analog" truncated Maxwellian flux distribution F (see section 1.3 above) with zero drift in π-direction but with different temperature [(T)\tilde]π ,([(T)\tilde]πT , [(V)\tilde] d = 0 ). This can be done by setting, e.g.,
 ~ v π = ~ α π · √ −ln(ξ) ·2 ,     ( note: ~ α π = ~ β π ,
 since
 ~ V d,π = 0 ),
(5.50)
ξ being a uniformly distributed random number between 0 and 1.
In order to find the weight correction factor associated with this procedure, we first note that the normalizing factor απ2 for the correct vπ distribution reads:
 απ2 = β2 ·exp(−M2π ) g(Mπ ) ,
(5.51)
with Mπ = Vd/ (√{ 2 } ·β) and g(x) = 1 + √{π} x (1 + erf(x)) ·exp(x2).
The conditions
 ~ M π = 0 , ~ α π = ~ β
then give the following result for the weight correction factor w ( [(v)\tilde]π ) to be applied for each sampled non-analog velocity component [(v)\tilde]π:
w (
~
v

π
) =
 Γ0Sh,ion ( ~ v π )

 ~ Γ 0Sh,ion ( ~ v π )
 =
 ~ β 2

β2
exp

 ~ v π Vd, π

β2
 ~ v 2π

2 β2

1 − β2

[(β2)\tilde]

· 1

g(Mπ)
.
(5.52)
This non-analog sampling of the source term introduces additional statistical fluctuations in the results, since the initial test particle weights fluctuate with variance
 σ2 = ⌠⌡ ∞ 0 d ~ v π w2 ( ~ v π ) · ~ Γ 0Sh,ion ( ~ v π ) − 1
 =
 g ⎛⎝ 2Mπ √b ⎞⎠ /g(Mπ)2 · 1 (2−b) b − 1
(5.53)
where b = 2 − β2 / [(β2)\tilde]. We use the ansatz b = C ·Mπ + 1 and compute numerically that value of C which minimizes σ2 in the range 0 ≤ Mπ ≤ 1. We find C = 0.6026. This means that the ratio of the "real" temperature T and the non-analog temperature [(T)\tilde] should be chosen to be T / [(T)\tilde] = 1 − 0.6026 ·Mπ. This choice optimizes the statistical performance. In principle, of course, any other [(T)\tilde] could have been used without introducing bias into the algorithm.
The mean energy of ions generated by this sampling procedure is (note: Γ0 is the distribution already normalized to flux one):
 E (Mπ)
 =
 mi 2 · ⌠⌡ ∞ −∞ ⌠⌡ ∞ −∞ ⌠⌡ ∞ 0 d3v′  v′2 ·Γ0Sh,ion (v′) =
 =
 mi 2 · ⌠⌡ ∞ −∞ dv1 ⌠⌡ ∞ −∞ dv2 ⌠⌡ ∞ 0 dvπ   (v12+v22+vπ2) ·Γ0Sh,ion (v′)
(5.54)
 =
 T ·γE ,
with the half sided integration for the π-component of v′, and:
 γE = ( Mσ2 + 2 ) + ⎡⎣ Mπ2 + 0.5· g(Mπ ) − 1 g(Mπ ) ⎤⎦ ,
(5.55)
We have decomposed the drift velocity Vd into a normal (to the surface) π-component Vd and a σ-component Vd parallel to the wall surface.
 Mσ = | Vd − Vd,π | √2 ·β = | Vd,σ | √2 ·β .
(5.56)
Note that for Mπ = 0 this reduces to the expected formula γE = 2 + Mσ2. Furthermore, if Mσ = 0 (normal incidence) and Mπ → 1 (Bohm sheath condition for Te = Ti), then γE for our forward sided fluxes approaches the value for total Maxwellian heat fluxes: γE = 2.449 + Mπ2 ≅ [5/2] +Mπ2, as it must.
We conclude that the above scheme is both an accurate and efficient random number generator for truncated shifted Maxwellian flux distributions of particles incident onto or emitted from surfaces. Of course, after having generated [(v)\tilde] with the weight w ([(v)\tilde] π) and before reflecting the particle via the kernel Cw , the acceleration of an ion in a sheath potential in direction π (normal to the target surface), as formalized in the sheath transmission kernel TSh mentioned above in Eq. (5.43), can easily be included.

1.5.1  The electrostatic sheath

The electrostatic sheath potential in front of a target surface is derived from assuming a given net electrical current jpl and a known ion particle current Γ+ (i.e., j+ = qiΓ+ electrical ion current, with qi=eZi being the ion charge, e the elementary charge and Zi the ion charge number) flowing over the target sheath edge. Both these currents, and the electron current j=qeΓ = −eΓ, are taken to flow parallel to the B-field at the plasma-sheath edge.
Often, but not necessarily, it is assumed: jpl = j+ + j = 0, i.e., locally ambipolar flow conditions.
We start with a Maxwell-Boltzmann density distribution ne(x) for the electrons, with electron temperature Te, in the electron retarding electric field with potential Φ(x). Here x is a coordinate parallel to the magnetic field B. The electron current density along the field is then (after reduction, possibly, due to secondary electron emission):
jT = −e ·1/4 ·ne,T ·veav ·(1−γs.e.e.) = −e ·1/4 ·ne,S ·exp[e∆ΦS,T/kTe] ·   ⎛

 8 kTe πme

·(1−γs.e.e.)
Here γs.e.e. is the secondary electron emission coefficient and veav is the average electron velocity. The subscripts S and T stand for the plasma-sheath interface and for the target surface, respectively and ∆ΦS,T is the (typically negative) voltage drop from the plasma sheath edge to the target.
In order to determine the sheath voltage drop ∆ΦS,T we equate the component normal to the target of this current density with the corresponding component of the net ion current density:
 −j−T ·cos(Ψ) = (j+T − jpl) ·cos(Ψ)
Hence the cosine of the angle Ψ of field line inclination against the target surface cancels out, as does the elementary charge e. The ion electrical current j+T ( = j+S, no ion sources within the sheath region) is given as
 j+T = e ∑ i Zi ·ni,S ·V||,i,S = e ∑ i Zi ·ni,S ·Mi,S ·ci,S.
The sum is over all ion species, Zi is the charge state of species i, ni,S is the corresponding ion density and V||,i,S the parallel (to the B-field) ion flow velocity of ion species i. Mi is the Mach number, and ci is the ion acoustic speed for species i.
Hence, the sheath potential difference ∆Φ in terms of given jpl, Te, γs.e.e. and ions flows Γi, reads:
 − e ∆Φ/kTe
 =
ln

 √ 2π

(1−γs.e.e.)

i

Zi ni,S

ne,S
·V||,i,S

[(me)/(kTe)]

jpl

e ·ne,S ·

[ˉ([(kTe)/(me)])]

(5.57)
 =
ln

 √ 2π

(1−γs.e.e.)

i

Zi ni,S

ne,S
·Mi,S ·ci,S

[(me)/(kTe)]

jpl

e ·ne,S ·

[ˉ([(kTe)/(me)])]

(5.58)
This expression (5.54) is evaluated by the function SHEATH of the EIRENE code, taking local plasma data for Te and ni and V||,i (ne from quasi-neutrality) at the point of impact of an ion at a wall boundary (e.g., the divertor target plate).
Special forms of (5.55) are often quoted in the literature. Writing for the ion acoustic speed ci explicitly
ci =   ⎛

 ⎛⎝ γkTi + kTe mi ⎞⎠

(γ denoting the adiabatic coefficient) and replacing this in (5.55):
ci ·   ⎛

 me kTe

=

γkTi

kTe
+ 1
· me

mi

1/2

we recover well known expressions for simpler cases. Consider, e.g., a single ion fluid case, ne = ni, with Zi = 1, and Mi = 1 (Bohm sheath condition for isothermal flow, γ = 1), furthermore zero net electrical current: jpl = 0 and zero secondary electron emission γs.e.e = 0. Then we find:
 − e ∆Φ/kTe = 1 2 ln ⎛⎝ 2πme mi ⎛⎝ Ti Te + 1 ⎞⎠ ⎞⎠
Hence, e.g., for a pure hydrogen plasma, and for TeTi, we have −e ∆Φ/kTe ≈ 2.5 and ≈ 2.84 for a pure deuteron plasma.

1.6  Combinatorial description of geometry

EIRENE is a Monte Carlo Code. Defining geometry in a Monte Carlo code is complex and laborious.
Of the three typical options for definition of the target geometry in Monte Carlo Codes:
• "Constructive solid geometry (CSG)",
• "boundary representation (BREP)"
• "Voxel Geometry"
EIRENE uses a combination of the latter two concepts.
In CSG a solid object is "constructed" from the intersection, union or other boolean operations of several "half-spaces". Each half space in specified by a relation F(x,y,z) ≤ 0, typically F is a low order polynomial. For example the neutronics code MCNP used e.g. for ITER nuclear safety studies, is based on CSG.
In BREP geometry, surfaces and volumes are build up from constituent parts. These parts may be related with algebraic functions (elementary BREPs) or complex functions such as Bezier functions or B-splines (advanced BREPS). Most CAD software uses BREPs in their native definition of geometry.
In the EIRENE BREP options surfaces are build from (bounded) algebraic functions. These surfaces are referred to as "additional surfaces", see section 2.3.2. EIRENE uses only algebraic functions, up to second order (some extensions to raise the order to fourth order have been started recently but are not ready to use).
In Voxel geometry, a 2D or a 3D object is constructed from identically shaped volume elements ("voxels"). In EIRENE these may be triangles or quadrangles in 2D, or tetrahedrons in 3D. New options to accommodate trilinear hexaedra as 3D voxels are currently being developed in connection with 3D edge codes such as EMC3. The bounding surfaces of these voxels are referred to as "standard grid surfaces" in EIRENE, see section 2.2. These volume discretisation meshes (VOXELS) can be read into EIRENE from external files. For example the mesh in B2 CFD plasma fluid code is read into EIRENE directly from the CFD grid generator in coupled B2-EIRENE applications. Wall meshes are build from these by assigning a surface reflection model to some of the (voxel-) surfaces, whereas internal (grid-) surfaces remain transparent for flow of test-particles.
An attempt has been made in 2008 to use CAD files directly for 3D wall surfaces. This was done via the already existing CAD-MCNP interface. We had then tried to turn the CSG input for MCNP (see above) into a BREP form for EIRENE. In principle this works. But in practice the de-tour via the CSG concept leads to prohibitively complex BREP input. A direct CAD-BREP interface seems now much more efficient to us and more natural anyway, but this has not been started yet in EIRENE due to lack of manpower.

A complete set of surfaces in an EIRENE geometry may consist of a combination of both: standard surfaces (hence: defining standard grid cells, (voxels), and additional surfaces, defining a general BREP discretisation (additional cells). A more complicated situation arises when additional surfaces intersect standard grid voxels. If these surfaces are transparent, they can be used to conveniently "measure" fluxes at any position. But if they are reflecting, or absorbing, e.g. to define a complex shaped bounding surface in an otherwise simple grid, then some care is needed because the automatic cell volume evaluation, which is available in EIRENE for standard grid cells (voxels), may not correspond anymore to the cell volume accessible for test particles. Cell averaged quantities (divided by cell volume) may then not necessarily be scaled properly. For a possible work-around see NLERG option in input block 1, section 2.1.
From the above it follows than in an EIRENE model surface elements building up the wall can be defined as bounded algebraic surfaces, with triangles as special case, i.e. as "additional surfaces, BREP", or can be bounding surfaces of the voxels, if the "grid of voxels" extends up to the real wall, or a combination of both.
As can be seen from the flowcharts, all the geometrical calculations needed for tracing test flights and for estimating volume- and surface averaged tallies are compiled in a "geometry block", which can easily be exchanged. In most cases (e.g., all installations of EIRENE outside FZ-Jülich) we load the most complete, fully 3D block "GEO3D", which allows fully 3 dimensional spatial resolution. Lower dimensional (0D, 1D or 2D) are contained as special cases, e.g. by choosing only grid surfaces with proper symmetry.
There is, in general, not a very large CPU benefit from running EIRENE with simpler specialized and geometry optimized blocks such as GEO2D or GEO1D for 2D or 1D models, respectively, because the saving in CPU costs in most cases is only around 10 - 20 %, so these options have not been supported since the late eighties of the last century anymore.
The following operations have to be executed by the block:
Suppose the three dimensional computational volume is discretized by a mesh of n surfaces Si, i=1,2,...,n . Each of these surfaces is defined by a set of Λi equations (unbounded surfaces):
 fλi (x,y,z) = 0 ,         λ = 1,2,..., Λi ,
(6.59)
and, optional, for each unbounded surface λ a subset of Πλ inequalities (e.g. for turning the unbounded surfaces into finite segments):
 gλ, πi (x,y,z) ≤ 0 ,         π = 1,2,..., Πλ .
(6.60)
This means that for all unbounded surfaces Si only the points, which fulfill all corresponding inequalities (6.57) are regarded by EIRENE as valid parts of the finite surface element.
Assume a test particle at position r0 = (X0,Y0,Z0) moving with speed unit vector Ω0 = (VELX,VELY,VELZ). The block then has to provide the nearest intersection r1 of the straight line ("sample trajectory")
 G(t): r0 + t Ω0
(6.61) amongst all possible intersections with the surfaces Si. Moving the "particle" from r0 to r1 and then repeating the block until the total track length l as sampled from the transport kernel T(l) (2.9a) is reached, allows to generate histories with all information available along the track to evaluate the estimators mentioned in section 1.2.
In EIRENE f and g are general 2nd order (with linear surfaces as special case) algebraic equations of the form
 h(x,y,z) = a0 + a1 ·x + a2 ·y +a3 ·z
 +
 a4 ·x2 + a5 ·y2 +a6 ·z2 +
 +
 a7 ·xy + a8 ·xz +a9 ·yz
(6.62)
Hence in the geometrical block nothing else but second order equations
t2 + p ·t + q = 0
have to be solved accurately and as efficiently as possible. We note in passing that each of this surfaces may be transparent, purely absorbing, semi transparent or (partially) reflecting. Surface sources may be defined on each surface. These surfaces may be "invisible" from one side and transparent, absorbing or reflecting from the other and each surface can "operate switches" if intersected by a test flight, such as abandoning the evaluation of atomic or molecular reaction rates (entrance into vacuum regions) along the flight, abandoning evaluation of volume averaged responses (e.g., in some uninteresting regions) or abandoning geometrical calculations (for those surfaces, for which it can be known a priori that they will not intersect a flight path originating at a particular surface or in a particular cell).

1.7  time dependent mode of operation

The section is currently rewritten. Please refer to the J.Nucl.Mater. paper: D.Reiter et al., Vol 220-222 (1995), 987-992.
The following issues have to be considered:
The coupling is explicit, hence CFL-type conditions for the time-step have to be obeyed. Here this means: the time-step must be shorter than the time between re-ionization and subsequent neutralization (in the volume or at the targets), because the fate of a particle during its intermediate ionized phase is not included in the explicit Monte Carlo procedure.
The noise level for collision estimators increases with smaller time-steps. Hence: avoid all collision estimators in t-dependent mode.

1.8  non-linear effects: coupling to plasma fluid models

The section is currently rewritten. Please refer to the report Jül-2872, Feb 1994, by G.Maddison and D.Reiter for details (ref. [9]).
Figure 1.1: Structure: EIRENE code
Figure 1.2: Structure: EIRENE code
Figure 1.3: Structure: EIRENE code
Figure 1.4: Convergence of B2/EIRENE, measured by B2-residuals. Method: implicit coupling, method of "saturated residuals", ref.[9]

1.9  non-linear effects: neutral-neutral collisions

Neutral-neutral collisions make the Boltzmann collision integral non-linear (bi-linear). More generally, inclusion of any binary collisions between two test-particle species i0, j0, elastic or not, have this effect. Important examples are, of course, elastic self collisions (i0,i0)el, e.g., H2H2 collisions leading to viscous effects in the gas, elastic cross collisions (i0,j0)el, such as HH2 collisions, often leading to a heating of the molecular gas due to hot charge exchange atoms. Furthermore there may be relevant inelastic binary particle collisions (i0,j0)inel, such as (H2H+2H+H3+) which is an important channel in hydrogen plasma chemistry, and (hν+ H(1s) → H(n=2)), the absorption of resonance Lyman-alpha photons under non-LTE conditions. In large dense divertor plasmas the ratio H(1s) to H(n=2) may have to be calculated consistently with the radiation field (photon Boltzmann equation), competing electron impact collisions and neutral atom transport.
Two methods of solution for this non-linear problem have been tested in EIRENE, 1) a parametrization of the single particle distribution functions involved in non-linear collision terms and an iterative scheme (see next subsection: BGK model collision terms), and 2) a Direct Monte Carlo Simulation (DMCS) technique based upon a swarm simulation, an operator splitting ("Streaming" and "Collisions") and a random simulation of binary collisions in the swarm. Currently only the first of these two options is still being used because, despite its approximate character, it proved to be much more effective than the full DMCS method, for the cases relevant to fusion research applications studied so far. A simple estimate for the relevance of neutral-neutral gas collisions is given by the so called Knudsen number Kn:
 Kn = λel l = 1 √2 σel nN l
(9.63)
Here λel is the mean free path for elastic neutral neutral collisions, σel is the typical cross section, nN the gas density (cm3) and l a characteristic dimension (cm). Fixing the hard sphere molecular diameter at d=4 ·10−8 cm results in a typical cross section of 1−2 ·10−15 cm2 to be used in this Knudsen number. For Kn << 1 a neutral fluid model (Euler, or Navier Stokes) is appropriate, while Kn >> 1 leads to the free molecular flow regime, in which only the linear collisions terms (collisions of neutrals with a given plasma host medium) are retained.

1.9.1  BGK approximation

See flowchart 9.1 (figure 1.12).
Further details can be found if ref. [18].

Consider first the elastic self collisions amonst particle species A. The BGK approximation to the Boltzmann collision integral reads:
 ∂fA ∂t =
(9.64)

1.10  Radiation transport, photon gas simulations

The theory of light and its interaction with matter is one of the best that physics can offer, and can predict virtually every observed phenomenon with great accuracy [Feynman, 1985, QED: The strange Theory of Light and Matter, Princeton University Press, Princeton]. In our simulations of radiation transport with the EIRENE code we assume a geometric optics model (i.e., essentially, "particle theory of light"), which is sufficient for all our purposes. Light is emitted, scattered, reflected and absorbed, but moves on straight lines between these events. I.e., a continuously varying index of refraction is not allowed, and we ignore most properties of light that depend on a wave or quantum model for their explanation (e.g. diffraction, or interference). However, the linear wave-optics phenomenon of polarization and the linear quantum optics phenomenon of fluorescence [absorption of a certain wavelength (line), and re-emission at another wavelength (line)] can be added to our linear geometric optics model.
We are concerned with the kinetic theory of interacting particles and photons, where the behaviour of each species is governed by a kinetic equation: the kinetic equation for photons being nothing else but the well known equation of radiation transfer, i.e., the "strangely normalized photon-Boltzmann equation". In contrast to plasma radiation (e.g., Bremsstrahlung) the radiation processes studied here are due to single, uncorrelated particles such that emission and absorption of the plasma is obtained by simple summation of all contributions of all individual particles. We therefore consider all particles as uncorrelated and in well defined one-particle states.
For continuum radiation (bound-free or free-free) the frequency of the photon in the atoms rest frame may be replaced by the frequency in the lab. frame, because all cross sections are slowly varying functions of frequency. By contrast, for line-radiation due to bound-bound transitions, the Doppler effect must be taken into account in all quantities that vary rapidly as function of frequency, such as cross sections, line-profile, or redistribution functions
The key atomic parameters for radiation transport are the Einstein (or Einstein-Milne) coefficients Amn, Bmn and Bnm of spontaneous emission, induced (stimulated) emission and absorption, respectively. We will discuss them here first for a gas of stationary atoms and include Doppler- and other line broadening effects later.
For a thermal radiation field (temperature T) in a gas of atoms with density n1 (lower state) and n2 (upper state) the rate equation expressing balance between emission and absorption reads:
 n1 B12 Rν(T) = n2[A21+B21 Rν(T)]
.
Rν(T) is a specific function characterizing the wavelength and temperature dependence of the thermal radiation field (i.e., the Planck-function). For example the specific intensity Iν (energy/area/time/solid-angle/frequency-interval) is often used, or the (angle)-integrated intensity Jν = ∫Iν dΩ, or the specific energy density uν = 1/c Jν (energy/volume/frequency-interval), etc. Also the wavelength dependence may be given either using frequency ν (Hz), ω (rad/s), wavelength λ or energy E = h ν = (h/2p) ω.
Depending upon all these choices the numerical factors in the definitions of the Bnm coefficients (sometimes even in the Anm) differ in the literature. For our purpose (Monte Carlo solution of the full radiation transfer equation) we need to define monocromatic mean free pathes for photons, rates and rate coefficients and cross sections for the individual processes.
In order to take advantage of the analogy between radiation and particle transport equations as much as possible, we will choose the numerical constants in the Einstein coefficients such that these parameters have the same meanings and even the same units in both cases of radiation and particle transport. I.e., we will use the "absorption coefficient" (better: "extinction coefficient", if scattering is allowed) αν = n σ ≡ Σt, (invers monocromatic mean free path), which directly corresponds to the "total macroscopic cross section" (dimension: 1/length) defined earlier for particle transport problems.

1.10.1  Line shape options

Emission of photons from a particular line is simulated by sampling firstly the location of emission, and secondly the frequency and direction of emission. Spatial distribution of photon emission is treated in EIRENE in exactly the same way as for any other particle type (neutrals, ions), and actually using the same modules of code for that. Emission from the transition ik from an upper state i=nl of an atom is treated, by abuse of language, as "volume recombination source": The position is sampled from the spatial distribution AikPnl([(r)\vec]), which, because the Einstein emission coefficient Aik is constant, is the same as the spatial distribution Pnl([(r)\vec]) of upper state particle density. ..... Given the birth point [(r)\vec] of the photon, and hence all local parameters of the host medium at this point, the frequency ν (in fact: the energy E=hν) is sampled from the normalized line shape function Φ(ν), .... The following line shape options are currently available:
• ISHAPE=0 δ-profile: mono-energetic emission (not yet available for absorption)
• ISHAPE=1 pure Doppler profile
• ISHAPE=2 Lorentz profile, truncated at zero
• ISHAPE=3 Doppler convoluted Lorentz profile, i.e. Voight profile
• ISHAPE=4 LORVDW: convoluted Lorentz and van der Waals pressure broadened profile
• ISHAPE=5 Doppler broadened LORVDW profile
• ISHAPE=6 normal Zeeman triplet, δ-profile in each component
• ISHAPE=7 Doppler broadened Zeeman-δ-profile in each component
• ISHAPE=8 Doppler broadened Zeeman-Stark-profile, 10 Lorentzian model, Rosato 2006

1.11  Charged Particle Transport

EIRENE also follows charged particles since EIRENE1987. The particle tracing is carried out in Subroutine FOLION. The differences of this module to the module FOLNEUT for neutral particles (and photons) are:
• Motion is not along straight lines between events, but instead determined by the electro-magnetic fields.
• Time stepping is not of Poisson type, but with constant increments, because the equation of motion has to be solved numerically.
• after each time-step a Coulomb collision kernel is modelled, i.e. a diffusive step in velocity space.
• after each time-step a diffusion-step is modelled, i.e. a diffusive step in physical space (not yet written).
• ionized test particles which hit a surface see an electrostatic sheath potential, by which they are either reflected, slowed down or accelerated, depending on the sign of their charge relative to the sheath potential.

1.11.1  Orbit integration

Particles following the B-field

A: no electric field

1.11.2  Coulomb Collision Models

Simple Coulomb Collision Models

Strongly simplified Coulomb collision models are usually sufficient for molecular ions, because of their very short lifetime until further fragmentation (dissociative excitation, dissociative recombination, etc...) to neutral products. In order to roughly account for thermalization effects and background plasma cooling or momentum transfer the energy, parallel and perpendicular velocity of such ions can be modified "continuously" during particle orbit integration, hence approximating the velocity space diffusion (Fokker Planck equation) by mean values of energy and momentum of the test particle. There are two classes of simple collision models: deterministic and stochastic.
A: Deterministic relaxation in velocity space   The energy and/or velocity of charged test particles are modified in a fully deterministic way along the particle trajectory. I.e.: along with orbit integration also energy, or certain velocity components are evolved according to ordinary (non-stochastic) differential equations, using appropriate relaxation times.
A1: Only energy relaxation
The simplest Coulomb Collision model is an averaged BGK-type relaxation of energy of test particles to the mean energy of the Maxwellian background ions. In the earliest versions of the EIRENE code (mid eighties) it was originally based upon the hydrocarbon ion equilibration model given by Langer, [20] and it is a slight generalization thereof. The equation for the energy εa of a test-ion a in a background of ions b with mean energy εb = 3/2 Tb is given as:
 dεa dt =−νεa,b ·εa     ,     εa(t=0)=εa0
(11.65)
This is a non-linear equation because the energy exchange collision rate (energy transfer rate) νεa,b depends on εa itself. The full expression for νεa,b is given below, it involves the Maxwell integral Ψ and its derivative Ψ′, see below.

In the limit of small test particle velocities (compared to thermal velocities of field particles b), i.e.: va << √{(2Tb/mb)}, i.e. xa/b = μbva2/2kTb << 1, as often appropriate for molecular ions, one has [21]:
 νεa,b = 2 νsa,b − ν⊥a,b − ν||a,b.
(11.66)
Here the energy transfer rate νεa,b is given via the "slowing down rate" νsa,b, the transverse diffusion rate νb and the parallel diffusion rate ν||b. In this limit on has for the basic "Spitzer rates" (loc. cit.):
 νsa,b
 =
 nb Za2 Zb2λa,b·6.8 ×10−8 μb1/2 μa ⎛⎝ 1+ μb μa ⎞⎠ Tb−3/2
(11.67)
 ν⊥a,b
 =
 nb Za2 Zb2λa,b·1.4 ×10−7 μb1/2 μa Tb−1/2ε−1
(11.68)
 ν||a,b
 =
 nb Za2 Zb2λa,b·6.8 ×10−8 μb1/2 μa Tb−1/2ε−1
(11.69)
and hence, using Equation (11.62) :
 νεa,b
 =
 nb Za2 Zb2λa,b·6.8 ×10−8 μb1/2 μa 1 Tb1/2 ⎡⎣ 2 1 Tb ⎛⎝ 1+ μb μa ⎞⎠ −2 1 ε a − 1 ε a ⎤⎦
(11.70)
 =
 νcε +νvε
(11.71)
with the "constant" and "variable" contribution νcε and νvε, respectively, defined as
 νcε
 =
 2 nb Za2 Zb2λa,b·6.8 ×10−8 μb1/2 μa 1 Tb3/2 ⎡⎣ 1+ μb μa ⎤⎦ = ~ ν ε ⎡⎣ 1+ μb μa ⎤⎦
(11.72)
 νvε
 =
 2 nb Za2 Zb2λa,b·6.8 ×10−8 μb1/2 μa 1 Tb3/2 ⎡⎣ − 3Tb/2 εa ⎤⎦ = ~ ν ε ⎡⎣ − 3Tb/2 εa ⎤⎦ .
(11.73)
Using this decomposition and definition of [(ν)\tilde]ε in equation (11.61) we find:
 dεa dt =− ~ ν ε ⎡⎣ 1+ μb μa ⎤⎦ ·εa + ~ ν ε · 3Tb 2
(11.74)
The solution of this linear differential equation is
 εa(t)=εa0 ·exp(− ~ ν ε [1+μb/μa] t)+ 3Tb 2 ⎡⎣ 1−exp(− ~ ν ε t) ⎤⎦
(11.75)
With μb=1 for protons, μa=16 for CH4, λa,b ≈ 10, and [1+μba] ≈ 1, hence [(ν)\tilde]ε ≈ 8.8 ×10−8 npTp−3/2 this is also the solution Eq. (11) given in [20] for thermalization of hydrocarbons in a Maxwellian hydrogen plasma. And this is exactly the approximation which was also used in EIRENE since about 1987 for all test ions traveling in a bath of background ions.
To be done: in case of drifting Maxwellian background: first transform into frame such that stationary Maxw. background, then relax test particle energy in this frame, then transform back...

The full expression for the energy relaxation frequency νεa,b, without the approximation of slow particles a, i.e. for all values of xa/b is:
 νεa,b(xa/b) = 2 [μa/μb Ψ(xa/b) − Ψ′(xa/b)] ν0a/b
(11.76)
where xa/b has been defined above, ν0a/b = ... and
 Ψ(x) =
2

 √ π

x

0
dy eyy
= erf(x) 2

 √ π
exx
(11.77)
 Ψ′(x) =
 dΨ/dx
= 2

 √ π
exx
(11.78)
With this general expression for the collision frequency Equation 11.61 can then only be integrated numerically along the particle trajectory.

1.11.3  The full Fokker-Planck Collision Model, A

The full Fokker-Planck Coulomb Collision Model, which was originally developed for a stand alone trace impurity ion transport code DORIS [22], is accommodated now as collision kernel in EIRENE: Thesis J. Seebacher, Univ. Innsbruck, soon to come

1.12  EIRENE flow charts

Figure 1.5: Structure: EIRENE code
Figure 1.6: Flowchart: Program EIRENE
Figure 1.7: Flowchart 2: Subroutine INPUT
Figure 1.8: Flowchart 6: Subroutine MCARLO
Figure 1.9: Flowchart 6.3: Particle loop in subr. MCARLO
Figure 1.10: Flowchart 9.1: Iterative mode for neutral-neutral interactions

Chapter 2 Description of formatted input file

General Remarks

All input data from unit IUNIN are read in subroutine INPUT using formatted I/O. There are 14 "blocks", each one starting with a card

*** n. TEXT

where n is the number of the input block, n=1,2,...,14. Some blocks are subdivided into so called ßub-blocks", each one starting with a card

** TEXT .

Each input block is described in one of the following 14 sections 2.1 to 2.14:

2.1 Input data for operating mode
2.2 Input data for standard mesh
2.3A Input data for "Non-default Standard Surfaces"
2.3B Input data for Ädditional Surfaces"
2.4 Input data for species specification and atomic physics module
2.5 Input data for plasma background
2.6 Input data for surface interaction models
2.7 Input data for initial distribution of test particles
2.8 Additional data for some specific mesh zones
2.9 Data for statistics and non-analog methods
2.10 Data for additional volumetric and surface averaged tallies
2.11 Data for numerical and graphical output
2.12 Data for plasma diagnostic module DIAGNO
2.13 Data for nonlinear and time-dependent mode
2.14 Data for interfacing with external codes

This last block 14 may be read in from the user supplied part (with any format specified there) from this same unit IUNIN, (but equally well from any other unit specified there), e.g., from the subroutine INFCOP for interfacing with other codes.

Units
cgs units are used, with two exceptions:

All temperatures and energies are in eV.
All particle fluxes are in Ampere (i.e. #/s ·1.602210−19), even neutral particle fluxes.
E.g., to convert a gas flow rate of x sl/min (standard liters per minute) first divide x by 22.4 (then: mol/min) next multiply by Avogadro's no. NA = 6.0221 ·1023, (i.e., then particles per minute), then divide by 60 (then: particles per second), and finally multiply by the elementary charge 1.6022 ·10−19 (then Ampere, as used for fluxes in EIRENE).
All densities are in cm−3.
All velocities are in either cm/s or in (isothermal) Mach-number units.
All geometrical data are in cm.
Note: consequently, e.g., energy densities are in eV/cm3, and energy fluxes are given in eV ·Amp, i.e., in Watt.
Hence: In order to convert from energy density units into pressure units in N/m2 (i.e., Pascal: Pa) one must multiply by 1.6022 ·10−19 (eV to Nm) and multiply by 106. Check: 1 Pa gas at 273 K corresponds to about 2.7 ·1014 ×0.026  eV/cm3 ≈ 7 ·1012  eV/cm3
For example, energy densities as computed by EIRENE defaults tallies 4, 5, 6 (see table 5.2) must be converted into "temperature densities" (a factor 2/3, if one neglects the energy of the directed motion) and then be re-scaled to pressure units as described above.
For further scaling to other pressure units note, e.g.:
101.3 Pa = 760 mTorr
Further: momentum is in g cm/s, and a momentum flux or momentum rate of change is in g cm/s ·Amp (per cell) and the momentum source density in the EIRENE tallies for interfacing with plasma codes is then in g cm/s Amp /cm3.
Hence, conversion from g cm/s Amp into Pa ·m2 (momentum source per cell in plasma momentum balance equations in SI units) is done by firstly multiplying with 10−5 (then: kg ·m/s ·Amp) and then dividing by 1.6022 ·10−19 (hence: kg ·m/s2). After dividing by the cell volume (m3) the momentum source is in Pa/m, or N/m3 (force density).

Input format

The format statement numbers in subroutine INPUT start with 666... for reading and with 777... for writing. The following conventions are used:
Standard Fortran naming conventions, i.e., variable names starting with the letters I,J,...,N are Integer, all other variables are Real.
In addition to these standard convention the following rules are employed:
The format for a card containing only real variables: (6E12.4)
The format for a card containing only integer variables: (12I6)
The format for a card containing only logical variables: (12(5L1,1X))
The format for a card containing text: (A72)
The format for a card containing K (K < 6) integer variables first and up to 6-K real variables next: (K(I6,6x),(6-K)E12.4)
An input card starting with '*' is recognized as a comment line, i.e., it does not contain input data. These comment lines should not start with '**' or '***', because in this latter case EIRENE would assume that the input for the current block is completed. It would expect the first card of the following block (or '* comment...' lines belonging to this next block) as next card. Any data not read from the input file, before the following block is read, are set to default values.
An arbitrary number of such comment lines may be included in the input file, but each block of comment lines must always be put at the beginning of a block or sub-block. I.e., comment lines are only identified by EIRENE if preceded by a card starting with '**' or '***' indicating a new sub-block or block respectively.

Example (see section 2.7 below):
*** 7. Data for initial distribution of test particles
* this is a comment
* there are NSTRAI=2 strata in this run
* this is a comment
2
** Data for first stratum
* this is a surface recycling source of helium ions, IPLS=2
* this is a comment
FFFF
.
.
.
.
** Data for second stratum
* this is a hydrogen ion volume re-combination source, IPLS=1
FFFF
.
.
.
*** 8. Data for some specific zones
.
.
.
.
and so on.

In very few cases formats other than these are used. They are then explicitly mentioned in this manual and in case of doubt the user should check in subroutine INPUT.

2.1  Input data for operating mode

General Remarks

The variables in this block control some of the more general options in EIRENE such as overall running time, usage of dump files but also a few parameters depending on the simulation model (drifts included or not, etc.).

The Input Block

TXTRUN
*** 1. Data for operating mode
first card: integer flags controlling global options for the entire run.
NMACH NMODE NTCPU NFILE NITER0 NITER NTIME0 NTIME
next input line is only for EIRENE2002 and younger, and optional in these versions (defaults given below).
NOPTIM, NOPTM1, NGEOM_USR, NCOUP_INPUT, NSMSTRA, NSTORAM, NGSTAL, NRTAL, NREAC_ADD
next line (second mandatory line): logical flags for global options.
NLSCL NLTEST NLANA NLDRFT NLCRR NLERG NLIDENT NLONE LTSTV
next input lines are only for EIRENE2005 and younger, and optional in these versions. There can be any number of these cards starting with character string CFILE in the input file.
"CFILE" DBHANDLE DBFNAME

Meaning of the input flags for operating mode

TXTRUN
At least one line of text to identify the run on printout and plot files.
Each text card must start with *. The first of these text lines will appear on all plots, the full text will be echoed at the beginning of the printout file
NMACH
Code-number for the computer used
= 1
CRAY
= 2
IBM
= 3
FACOM
= 4
VAX
This flag is currently not in use. It only affects the round-off error margin used for geometrical computations. In case NMACH=1, the error tolerated in the geometry routines is EPSGEO=1.D-10, and EPSGEO=1.D-6 for all other values of NMACH.
NMODE
Operating mode
= 0
EIRENE run as stand-alone-code. Code coupling segment  couple_Dummy.f  may be used. Input block 14 has a fixed format, see section 2.14.
≠ 0
EIRENE calls the code coupling subroutine INFCOP in the code coupling segment  couple_Name.f  where 'Name' is a character string identifying the particular external code, to which EIRENE is to be coupled (e.g.: Name=B2, Name=DIVIMP, Name=U-file, Name=FIDAP, Name=EMC3, Name=OSM, etc.).
Hence: the routine INFCOP is called by EIRENE for communication with external data sources, e.g., other codes. The calling program for the first 3 entries of INFCOP is subroutine INPUT, and the call to subroutine INFCOP is after reading 13 blocks from the formatted input file (READ (IUNIN,...)). A 14th block of the input file is read from subroutine INFCOP with the format (if any) specified there. At entry IF0COP geometrical data are provided (overruling corresponding data in input block 2).
At entry IF1COP plasma profiles are defined (overruling the background data specified in block 5).
Any other data (e.g., surface- or volume source distributions overruling the data in input block 7) are expected from entry IF2COP.
> 0
If NMODE > 0, subroutine INFCOP is called once again after each completed stratum (at the entry IF3COP(ISTRAA,ISTRAE)) to return data to another code (post-processing). The call to IF3COP is from subroutine MCARLO, at the end of the DO 1000 -loop over the strata. One final call to the interfacing routine (entry IF4COP) is implemented after the calls to the EIRENE printout- and plot routines, e.g. to perform global balances etc. after summation of the contributions from the individual strata.
NTCPU
Maximum number of CPU seconds allowed for this run. NTCPU must be less than or equal to the time parameter in the job-card (if any).
If more than one iteration is carried out (NITER, see below, this section) or more than one time cycle (NTIME, see below, this section), then each iteration or cycle can take up to NTCPU seconds.
In 2008 the definition of NTCPU was slightly changed: the initialisation overhead is not any longer included in NTCPU. I.e., this variable NTCPU now is close to the true Monte Carlo sampling time. The total cpu-time, including initialisation overhead as well as post-processing is printed at the end of an EIRENE run (see: "total cpu-time of this run" in printout file).
NFILE
Flag for the use of dump files FT10, FT11, FT12, FT13, FT14 and FT15.
= 0: Neither reading nor saving of data is done.
= JKLMN:
N = 1
EIRENE writes output data from this run into files FT10 and FT11 to save them for other plot or printout options in later 'read only' runs with the same input file
N = 6
same as N=1, but only the data for the sum over all strata are written (sufficient, e.g., for BGK-iterations, or for post processing, (subroutine DIAGNO, see input block 12) as long as this post processing only is to be done on total results (not on contributions from individual strata).
N = 2
EIRENE reads output data from an earlier run from files FT10 and FT11. No new particle histories are computed, only the requested output is printed and plots are produced. Iteration cycles or time cycles (if any) are abandoned.
N = 7
same as N=2, but only the data for the sum over strata are read. See N=6 option described above.
M = 1
EIRENE writes geometry data into file FT12. These data are the output from subroutines GRID and VOLUME in the initialization phase.
M = 2
EIRENE reads geometry data from file FT12. The geometry subroutines GRID and VOLUME are not called. This option should be used if the geometry has not changed as compared to an earlier run. It reduces the CPU costs for the overhead.
L = 1
EIRENE writes plasma data and atomic and molecular data ("A&M data") into file FT13. These data are output from subroutines PLASMA, XSECTA, XSECTM, XSECTI, XSECTP in the initialization phase, or at the end of a run, if some data have been modified in subroutine MODUSR, see NITER flag below. Also the primary source parameters (input block 7) are saved, for later iterations or time-steps.
L = 2
EIRENE reads plasma data, A&M data from file FT13. It also reads the entire input block 7 from FT13, and overwrites the input read from this block 7 by those data. Routines PLASMA, XSECTA, XSECTM, XSECTI, XSECTP are not called. This option should be used if neither the plasma background nor the selection of atomic processes to be used, nor the primary source model (block 7) has changed as compared to an earlier run. It reduces the CPU costs for the overhead.
L = 3
Acts as if both NFILE-L=1 and NFILE-L=2. For continuation of iterative calculations, for example.
L = 4
Same as NFILE-L=3, except: primary source data (input block 7) are not read from fort.13, but are taken as specified in input block 7. This allows to modify these primary sources parameters during iterations or time-steps not only via module MODUSR, but directly via input file.
L = 6,7,8,9
Same as L=1,2,3,4, respectively, but XDR file format is used.
K = 1
EIRENE saves some data for optimizing non-analog sampling and stratified source sampling on file FT14.
K = 2
EIRENE reads from file FT14 and optimizes operation. Currently only allocation of CPU time in stratified source sampling is optimized.
K = 3
Acts as if both NFILE-K=1 and NFILE-K=2.
J = 1
Only for time-dependent option (see block 13). EIRENE writes "census data" at the end of last time-step onto file FT15.
J = 2
EIRENE reads "census data" from file FT15, and uses it as initial distribution for the coming next time-step.
J = 3
Acts as if both NFILE-J=1 and NFILE-J=2.
NITER0
Initial iteration number: (Default: NITER0=1) Irrelevant parameter. Only needed for book keeping and printout. EIRENE labels the iterations from NITER0 to NITER.
NITER
Number of iterations, if EIRENE runs in "iterative mode".
> 0
EIRENE calls user supplied subroutine MODUSR after completing the run; some model parameters may be modified here for the next iteration step, and some results from the previous step may be saved on a file.
> 1
EIRENE recalls itself but does not read from the formatted input file again. This recalling is repeated NITER times (including the first iteration). The CPU time NTCPU is used for each iteration. Hence the true CPU time then is NTCPU ·NITER .
NTIME0
Initial time-cycle number: (default: NTIME0=1). Irrelevant. Only needed for book keeping and printout. EIRENE labels the time steps from NTIME0 to NTIME.
NTIME
Total number of iterations ("cycles") in time carried out in one single run. The total time per cycle is defined as NTMSTP · DTIMV (see below, input block 13).
After each time-cycle the subroutine TMSTEP is called. In this routine the "census arrays", which store the test particle population at time ti:
ti = ti−1 + ∆t = t0 + it = TIME0 + ITIME·[NTMSTP ·DTIMV]
are filled and prepared for the next time-cycle. The census arrays from the previous time cycle (if any), i.e., at t = ti−1, are overwritten here.
After the last time-cycle, the census arrays are written on file fort.15, in order to permit continuation in time in a next run.
The background conditions (and source distributions or any other input parameters) for the next time cycle can be modified in subroutine TMSUSR(ti), which is called from subroutine TMSTEP.
If the population on the census array is not empty or known from a previous run (NFILE-J flag, see above), then this census population defines one additional stratum for the current cycle. I.e., the census population then determines the initial condition for the distribution function f(r,v,i,t=t0). The source strength FLUX is computed from that initial condition.
NOPTIM
Default: 1
NOPTM1
Default: 1
NGEOM_USR
for value =1: user defined (external) geometry package (LEVGEO=10), then no grid storage is provided in EIRENE. Default: 0
NCOUP_INPUT
=1: Storage for data transfer via coupling routines, =0: no such storage. Default: 1
NSMSTRA
Sum over strata disabled for value 0, enabled for value 1. Default: 1
NSTORAM
Storage vs. speed in atomic data evaluation. Maximum storage, fastest computation: =9. Minimum storage, maximum work (slowest option) =0. Default: 9
NGSTAL
Storage for spatial distribution of surface tallies on non-default standard surfaces for value =1. For value =0: only total (spatially integrated) surfaces fluxes. Default: 0
NRTAL
Condensing mesh cells into fewer larger cells, for output volume tallies. The underlying fine mesh has NRAD cells (see input block 2). The coarser mesh, obtained from condensing cells into one larger cell, has NRTAL cells. Default: 0: Then internally: NRTAL=NRAD, and no condensation is carried out.
NCLTAL(IRAD) = IRTAL: grid cell IRAD is condensed into the larger cell IRTAL. I.e.: output is average over a larger cell IRTAL, which is comprised of all cells IRAD such that NCLTAL(IRAD) = IRTAL, IRAD=1, NRAD
The input data in input block 5 (background medium) are always given on the fine mesh (size: NRAD)
For output tallies (cell averaging) several cells IRAD1, IRAD2,... can be condensed into one larger cell IRTAL. Cell volumes, scoring, statistics are automatically done on the coarser grid (size: NRTAL).
The index array NCLTAL(NRAD) can be defined in the problem specific "user" routines (see section 3), e.g.: INIUSR, GEOUSR, etc..
Storage for additional reaction decks read onto EIRENE arrays in USR-routines, post-processing, etc. Default: 0
NLSCL
Some volume averaged tallies are re-scaled in order to exactly preserve the total number of particles, which otherwise would be the case only up to statistical precision (due to the use of track-length estimators). EIRENE computes three factors FATM, FMOL and FION such that particle balances for atoms, molecules and test ions, respectively, are accurately observed, if NLSCL = TRUE.
NLTEST
Tests for consistency between cell numbers and geometrical data along the particle tracks are carried out at each point of collision. If inconsistencies are detected, the history is stopped and an error message is printed. The contribution of these particles to the particle- and energy balances is stored in the bins "PTRASH" and "ETRASH" respectively.
NLANA
De-activates (NLANA=.TRUE.) all non-analog sampling distributions, such as biased source sampling, splitting, etc.. Select NLANA=.TRUE., if particle trajectory plots are used to get an intuitive picture of what is going on physically.
NLDRFT
Drift component in the bulk ions velocity distribution is included. (see input block 5). Otherwise an isotropic Maxwellian distribution is assumed for the background particles.
NLCRR
Correlation sampling is used.
NLERG
The case is automatically reduced to a case for estimating cell volumes by utilizing an "ergodic property".
NLIDENT
In multi-processor calculation mode NLIDENT forces all processors working on the same stratum to use the same sequence of random numbers and hence to carry out exactly identical work. This flag can be used to test the parallelized code version.
NLONE
The case is automatically reduced to a single species and "one speed transport" problem. (not ready, don't use)
LTSTV
free flag used for testing purposes
Optional input cards, for pathways to external databases. Their presence is identified by the starting character string CFILE of these cards.
DBHANDLE
Name of the database which is specified in this "CFILE" card. The format of a data file is identified in EIRENE by its name. Possible database names currently recognized are:
AMJUEL, HYDHEL, METHAN, H2VIBR, HYDRTC, ADAS, SPECTR, PHOTON, POLARI,

SPUTER, TRIM, TRMSGL,
with volumetric collision databases in the first block, and surface interaction databases in the second.
DBFNAME
Absolute or relative path and name of file to be used for this database.
Example:
CFILE AMJUEL /home/path/AMdata/amjuel.tex
CFILE TRIM /home/path/Surfacedata/TRIM/trim.dat
Default: If the specification of an external database name and path, which is later referred to in the run, is not given here, then these data files are expected, usually with the same name, in the directory in which the calculation takes place. Due to historical reasons and for backward compatibility for some data files also different default names are recognized.
Default file names to be used in this directory currently are
AMJUEL,
HYDHEL,
METHAN or "METHANE", for database name METHAN
H2VIBR,
HYDRTC,
SPECTR,
PHOTON,
POLARI,
SPUTER,
TRIM or "fort.21", for database name TRIM
For ADAS (AMdata) and TRMSGL (Surfacedata) there is an exception from this scheme. In these cases no default file names are available. And input variable DBFNAME specifies only the path to the directory which is the root of the ADAS and TRIM data tree, respectively.
Example:
CFILE TRMSGL /home/path/Surfacedata/TRIM/
The particular data file from these trees (e.g.: SCD96 in case of ADAS, or A_on_B in case of TRMSGL), the chemical element and the charge stage of a particular ion Z=1,... Z+1 (e.g.: C, 4), are detailed later in input block 4 (for volumetric collision processes based on ADAS adf11 format), and input block 6 (TRMSGL, for the surface reflection database based on TRIM) see below, sections 2.4, 2.6, respectively.

2.1.1  The NLERG option for cell volumes

All non-transparent surfaces are made perfectly reflecting (mono-energetic, cosine distribution). All volumetric collision processes are de-activated. Except, possibly, for the very first flight, there is only either one atomic species (IATM=1) or one molecule species (IMOL=1). Unless otherwise specified in input block 13, a single time-step, with time limit DTIMV=1.0 (s) is set and up to NPRNLI=10 test flights are launched. Under such conditions, the particle density should be constant throughout the grid. Statistically significant deviations from that constancy indicate wrong cell volumes as compared to the volumes seen by the test particles (e.g. due to additional surfaces intersecting the grids, or complicated shapes of some additional cells (see input blocks 2 and 3). In a second run the volumes of these cells can then be set explicitly, e.g. in input block 8.

2.1.2  The NLMOVIE option for making movies of trajectories

documentation to be written

2.2  Input for Standard Mesh

General Remarks

EIRENE follows atoms, molecules or test ions in a 3-dimensional computational box. This box is discretized by zones (mesh cells), the boundaries of which are defined either by regular "standard mesh surfaces", i.e., co-ordinate surfaces, which are described here, and/or by zones whose boundaries are defined more generally by "additional surfaces" (see below block 3B).
As long as test particles travel inside the "standard mesh", all cell indexing computations are done automatically. At transition into the more general "additional mesh" and inside these additional cells, this indexing has to be specified by the flags ILSWCH, ILCELL etc. by the user, in input block 3B.
There are up to 3 sets of standard surfaces defined by the arrays RSURF, PSURF and TSURF.
Here RSURF(1), RSURF(2), ..., RSURF(NR1ST) are grid-points in the radial or x direction (depending on the geometry level, see below) defined by the input data in sub-block 2A. I.e., there are NR1STM= NR1ST-1 cells in radial- or x direction. With regard to the terminology in 1.6, eq. 6.56, a corresponding co-ordinate surface is labelled with a double index 1,IR and is defined by the equation
f11,IR (x,y,z) =constant=RSURF(IR).
PSURF(1), PSURF(2), ..., PSURF(NP2ND) are grid-points in the poloidal or y direction (sub-block 2B), and there are NP2NDM=NP2ND-1 cells in this direction. A corresponding co-ordinate surface is labelled 2,IP and is defined by the equation
f12,IP (x,y,z) =constant=PSURF(IP).
In case of the NLPLG option ( (mesh in x-y plane generated by a set of polygons) for simplicity and by an abuse of language we sometimes refer to the polygons themselves as radial surfaces or "radial polygons" RSURF and refer to the polygons obtained by connecting all points with a particular index along the various radial polygons as "poloidal polygons" PSURF. No use is made in EIRENE of the arrays RSURF and PSURF in connection with the polygon mesh option. The NLPLG option permits to run EIRENE on 2D computer generated meshes with grid cuts in one direction.
Furthermore, TSURF(1), TSURF(2), ..., TSURF(NT3RD) are grid-points in toroidal or z direction (sub-block 2C). A co-ordinate surface is labeled 3,IT and is defined by the equation f13,IT (x,y,z) =constant=TSURF(IT).
In addition to these grids and, depending on the geometry options chosen, EIRENE then computes further "derived" grid arrays, for example a "(radial) surface labeling grid" RHOSRF, which may or may not be identical to the RSURF grid, and a "zone centered (radial) surface labeling grid" RHOZNE(IR) defined by:
RHOZNE(IR) = 0.5 · (RHOSRF(IR)+RHOSRF(IR+1))
These grids may be used for the input of closed form radial- or x direction plasma profiles or as abscissae in plots of e.g. poloidally (or y) and toroidally (or z) averaged tallies.
We have found it convenient to allow for several (NBMLT) identical copies of such "standard meshes", separated and bounded by arbitrary "additional surfaces". Input flags for this "mesh multiplication" option are comprised in sub-block 2D.
As mentioned above, in addition to these regular grids one can define cells of almost arbitrary complexity by using the "additional surfaces" defined in sub-block 3B. Such surfaces are labeled 0,IA in EIRENE. Appropriate cell number switching flags must be set on transparent additional surfaces separating these general "additional cells". Data relevant for the additional cells (e.g. the volume) are read in sub-block 2E. The number of additional cells is NRADD.
Therefore, the total number of cells in an EIRENE run is NSBOX, with
NSBOX = (NR1ST · NP2ND · NT3RD) · NBMLT + NRADD .
Cell volumes in the standard mesh region, as well as the center of mass vector in each of these cells are computed automatically. The corresponding data in the additional cell region must be specified by the user, e.g. in the input block 2E or 8, or in the user supplied routine GEOUSR (section 3.6), which is called by EIRENE in the initialization phase.

The Input Block

*** 2. Data for Standard Mesh
 (INDGRD(J), J=1,3)
** 2A. x- or radial co-ordinate surfaces
NLSLB  NLCRC  NLELL  NLTRI  NLPLG  NLFEM  NLTET NLGEN
NR1ST  NRSEP  NRPLG  NPPLG  NRKNOT
IF (INDGRD(1).LE.5) THEN
IF (NLSLB.OR.NLCRC.OR.NLELL.OR.NLTRI) THEN
RIA  RGA  RAA  RRA
IF (NLELL) THEN
EPIN  EPOT  EPCH  EXEP
ELIN  ELOT  ELCH  EXEL
IF (NLTRI) THEN
TRIN  TROT  TRCH  EXTR
ENDIF
ENDIF
ELSEIF (NLPLG) THEN
XPCOR  YPCOR  ZPCOR  PLREFL
(NPOINT(1,K) NPOINT(2,K), K=1,NPPLG)
DO 21 I=1,NR1ST
(XPOL(I,J) YPOL(I,J),  J=1,NRPLG)
21       CONTINUE
ELSEIF (NLFEM) THEN
XPCOR  YPCOR  ZPCOR
NRKNOT
(XTRIAN(I),I=1,NRKNOT)
(YTRIAN(I),I=1,NRKNOT)
DO ITRI=1,NR1ST
I,NVERT(1,ITRI),NVERT(2,ITRI),NVERT(3,ITRI),
NEIGH(1,ITRI),NSIDE(1,ITRI),IPROP(1,ITRI),
NEIGH(2,ITRI),NSIDE(2,ITRI),IPROP(2,ITRI),
NEIGH(3,ITRI),NSIDE(3,ITRI),IPROP(3,ITRI)
ENDDO
ENDIF
ENDIF
ENDIF
** 2B. y- or poloidal grid surfaces
NLPOL
NLPLY  NLPLA  NLPLP
NP2ND  NPSEP  NPPLA  NPPER
IF (INDGRD(2).LE.5) THEN
YIA  YGA  YAA  YYA
ENDIF
** 2C. z- or toroidal grid surfaces
NLTOR
NLTRZ  NLTRA  NLTRT
NT3RD  NTSEP  NTTRA  NTPER
IF (INDGRD(3).LE.5) THEN
ZIA  ZGA  ZAA  ZZA  ROA
ENDIF
** 2D. mesh multiplication
NLMLT
NBMLT
IF (NBMLT GT 1) (VOLCOR(NM), NM = 1,NBMLT)
** 2E. additional cells outside standard mesh



Meaning of the Input Variables

INDGRD
This index controls the meaning of input variables of different standard grid options.
INDGRD(1) for the radial (or x-direction) grid
INDGRD(2) for the poloidal (or y-direction) grid
INDGRD(3) for the toroidal (or z-direction) grid
Data for first standard mesh: radial or x grid RSURF

A radial or x grid is defined. Otherwise the complete sub-block 2A may be omitted. Depending upon the logical parameters in the next input card the "geometry - level" variable LEVGEO is set internally.
LEVGEO=1
cartesian coordinates x (and y)
LEVGEO=2
polar coordinates r (and θ)
LEVGEO=3
general curvilinear coordinates: a full 2D mesh (polygonal coordinate lines) is used in the x - y plane. Grid cuts are permitted in the y-direction.
LEVGEO=4
a 2D "finite element" mesh of triangles is used in the x - y plane
LEVGEO=5
a 3D "finite volume" mesh of tetrahedrons used
LEVGEO=6
a general, user defined geometry block is used. All geometrical calculations are performed in problem specific routine VOLUSR, TIMUSR, ...etc.
If NLRAD=.FALSE., then no spatial grid is defined and the default geometry level LEVGEO = 1 is used. Volume discretisation may still be achieved "by hand" by defining "additional surfaces" (input block 3b) and appropriate cell number switching.
INDGRD(1)
= 1
standard grid option; the input parameters are used as described below.
= 2
As for INDGRD(1)=1, but in case of LEVGEO = 2 by this option a radial grid is defined such that the spacing of radial surfaces is not equidistant at some poloidal position, but instead such that the area enclosed between two neighboring surfaces is kept constant. The input variable RGA is irrelevant in this case.
= 3,4,5
not in use
= 6
Data for NR1ST radial surfaces are read from the work array RWK. These data must have been written onto RWK in the user supplied subroutine INFCOP. By this option grid parameters may directly be transferred into EIRENE from other files, e.g. from plasma transport codes (see below: section 2.14. Data for interfacing routine INFCOP).
One and only one of the next 6 logical variables must be .TRUE.
NLSLB = .TRUE.

Geometry level: LEVGEO = 1
Cartesian geometry, the x co-ordinate is discretized by setting
xI=RSURF(I) I=1,NR1ST .
Furthermore, the flux-surface labeling grid RHOSRF(I) is identical with the grid RSURF(I).
NLCRC = .TRUE.

Geometry level: LEVGEO = 2
Cylindrical or toroidal geometry, mesh of concentric, circular surfaces. Polar coordinates are used in the x-y plane. The radial surfaces are given by r2 = x2 + y2 = const. and r is discretized by setting rI=RSURF(I) I=1,NR1ST. Furthermore RHOSRF(I) = √{ [(AREA)/(π)] } where AREA is the area inside surface number I. Thus for this option one has again: RHOSRF(I) = RSURF(I), I=1,NR1ST .
NLELL = .TRUE.

Geometry level: LEVGEO = 2
Mesh of nested, but not necessarily concentric or confocal elliptical flux surfaces. The equation for the "radial" surface is (xEP)2 + ( y/EL )2 = r2. r is discretized by setting rI=RSURF(I) I=1,NR1ST.
EP and EL may vary with r. These parameters are stored in the arrays EP(I),EL(I), I=1,NR1ST which now are used in addition to RSURF to define one surface.
RHOSRF: as in NLCRC option.
Note: RHOSRF and RSURF may differ in this case.
NLTRI = .TRUE.

Geometry level: LEVGEO = 2
to be written: triangularity in mesh of nested closed algebraic surfaces
NLPLG = .TRUE.

Geometry level: LEVGEO = 3
The mesh in the x-y plane is described by NR1ST polygonal arcs of length NRPLG each. A polygon may consist of several "valid" and "invalid" parts (to account for "grid cuts" in CFD meshes). The "invalid" parts of a polygon are not seen by test particles and are allowed for in EIRENE only in order to facilitate index mapping in case of runs coupled to plasma transport models, which resort to computer generated meshes including grid cuts.
The polygons must not intersect each other.
In this case RHOSRF(1)=0., and RHOSRF(I) is the area enclosed by polygon number 1 and polygon number I.
NLFEM = .TRUE.

Geometry level: LEVGEO = 4
The mesh in the x-y plane consists of NR1ST triangles, composed from NRKNOT knots.
In this case a flux surface labeling gird RHOSRF is not defined.
NLTET = .TRUE.

Geometry level: LEVGEO = 5
3D discretisation of volume by tetrahedrons. For this grid option please make contact to the authors.
NLGEN = .TRUE.

Geometry level: LEVGEO = 6
Arbitrary geometrical configuration. Mesh consists of NR1ST arbitrarily shapes cells (in any dimension). Particle tracing routines must be provided by user (VOLUSR, SAMUSR, TIMUSR, LEAUSR)
NR1ST
Number of grid-points in the radial (or x-direction) standard mesh
if NR1ST ≤ 1, no radial (or x-direction) standard mesh is defined.
NRSEP
This flag is active for LEVGEO = 1 or LEVGEO = 2. Otherwise it is irrelevant.
The standard mesh is composed by two equidistant x- or radial grids of co-ordinate surfaces with different grid density. There are NR1ST-NRSEP+1 grid-points in the first, and NRSEP grid-points in the second part. The grid-point RSURF(NR1ST-NRSEP+1) belongs to both parts.
RIA
left endpoint of standard grid (internally set ≥ 0 if LEVGEO =2); RSURF(1)=RIA
RGA
boundary separating first and second part of standard grid with different grid-point densities; RSURF(NR1ST-NRSEP+1)=RGA
RAA
right endpoint of standard grid: RSURF(NR1ST)=RAA
RRA
if RRA > RAA, one additional zone is defined
RSURF(NR1ST+1)=RRA, and NR1ST is increased by one.
(irrelevant, if RRA ≤ RAA)
if NLELL = .TRUE. :
EPIN
Value of EP(r) for cylindrical co-ordinate surface number 1 with r1=RIA (see: NLCRC = .TRUE. option above)
EPOT
Value of EP(r) for cylindrical surface number NR1ST with
rNR1ST=RAA
EPCH
Value of EP(r) for cylindrical surface number NR1ST+1 with
rNR1ST+1=RRA
(irrelevant if RRA ≤ RAA)
EXEP
The variation of the "shift function" EP(r) with r is given by
EP(r) = EPIN + ( [(rRIA)/(RAARIA)] ) EXEP·(EPOTEPIN)
ELIN
Value of EL(r) for cylindrical surface number 1 with
r1=RIA (see: NLELL = .TRUE. option)
ELOT
Value of EL(r) for cylindrical surface number NR1ST
with rNR1ST=RAA
ELCH
Value of EL(r) for cylindrical surface number NR1ST+1 with
rNR1ST+1=RRA
(irrelevant for RRA ≤ RAA)
EXEL
The variation of the "ellipticity function" EL(r) with r is given by
EL(r) = ELIN + ( [(rRIA)/(RAARIA)] ) EXEL·(ELOTELIN)
if NLTRI = .TRUE. :
to be written
if NLPLG = .TRUE. :
NR1ST
number of polygons for discretisation in "radial" or x direction
NRPLG
Number of points per polygon
NPPLG
Number of valid parts on each polygon. Each polygon is described by the x and y co-ordinates of NRPLG points. It is not necessary that all this points are used for the polygon. One can cut the polygon into several valid parts interrupted by parts which are not seen by the test particles. ( Default : NPPLG = 1 ). This option facilitates the use of 2-d computer generated meshes which contain topological grid cuts.
XPCOR,YPCOR

shift whole mesh by that vector in x,y-plane
RFPOL
if RFPOL > 0., one additional polygon zone is defined, at a distance RFPOL outside the polygon NR1ST, and then NR1ST is increased by one.
(irrelevant, if RFPOL ≤ 0.)
NPOINT(1,J)

Index of the first point of the valid part number J (same for each radial polygon)
( Default : NPOINT(1,1) = 1 )
NPOINT(2,J)

Index of the last point of the valid part number J (same for each radial polygon)
( Default : NPOINT(2,1) = NRPLG )
XPOL(K,I)
x-co-ordinate of the polygon point number K on polygon number I
YPOL(K,I)
y-co-ordinate of the polygon point number K on polygon number I
if NLFEM = .TRUE. :
NR1ST
number of triangles for discretisation in x-y plane
XPCOR,YPCOR

shift whole mesh by that vector in x,y-plane
NRKNOT
There are NRKNOT knots, by which the triangles are defined
XTRIAN,YTRIAN

x and y co-ordinates of the knots, respectively.
NVERT(I,ITRI)

Each triangle ITRI is defined by 3 points P1, P2 and P3 from the set of NRKNOT knots.
NVERT(I,ITRI) is the number of point PI (I=1,2,3) in the set of knots for triangle ITRI.
NEIGH(I,ITRI)

The three sides of each triangle S1, S2, S3 are defined such that S1 connects P1 and P2, S2 connects P2 and P3 and S3 connects P1 and P3.
NEIGH(I,ITRI) is the number of the neighboring triangle at side SI for I=1,2,3.
NSIDE(I,ITRI)

NSIDE(I,ITRI) is the number (1,2 or 3) of the side of the neighboring triangle, which corresponds to side SI of the triangle ITRI to be written
IPROP(I,ITRI)

ISTS = ABS(IPROP) is the integer, by which a particular surface property is assigned to side I of the triangle ITRI. ISTS=0 stands for default grid option (transparent surface, cell indexing is done automatically). Otherwise ISTS is the number of an additional (ISTS < NLIMI) or non-default standard (NLIM < ISTS < NLIM+NSTSI) surface option as read in sub-blocks 3a or 3b, respectively. By default the normal vector of each side of a triangle points out of the triangle. In case IPROP < 0, this vector points into the triangle. This is relevant at surface options with ILSIDE ≠ 0.
if NLTET = .TRUE. :
Documentation not available here. Make contact with the authors.
if NLGEN = .TRUE. :
NR1ST
number of cells in otherwise arbitrary mesh option NLGEN.
Data for second standard mesh: poloidal or y grid PSURF
NLPOL
A poloidal or y grid is defined. Otherwise the complete block 2B may be omitted and the volume averaged tallies are then automatically integrated over this co-ordinate.
INDGRD(2) = 1

standard grid option
INDGRD(2) = 2,3,4,5,6

not in use (default: INDGRD(2)=1)
NP2ND
Number of grid-points in y- or poloidal direction
YIA,YGA,YAA,YYA

for the LEVGEO=1 option: the y grid PSURF is defined in the same way as the x grid, using the parameters YIA,YGA,.... (cm) instead of RIA,RGA,....
for the LEVGEO=2 option: the θ grid PSURF is defined in the same way as the r grid, using the parameters YIA,YGA,...(poloidal angle in degree) instead of RIA,RGA,...
for all options LEVGEO > 2: this input card is irrelevant.
Defaults: (needed for scaling, i.e., cell volumes):
YIA=0., YGA=0., YAA=1., YYA=1. in case of LEVGEO=1,
and YIA=0., YGA=0., YAA=360., YYA=360. in case of LEVGEO=2.
Data for third standard mesh: toroidal or z grid TSURF
NLTOR
A toroidal or z grid is defined. Otherwise the complete block 2C may be omitted and the volume averaged tallies are then automatically integrated over this co-ordinate.
In case NLTOR = TRUE, sub-block 2C must be read.
INDGRD(3) = 1

standard grid option
INDGRD(3) = 2,3,4,5,6

not in use (default: INDGRD(3)=1)
One and only one of the next four boolean variables must be TRUE
NLTRZ=TRUE

cylindrical approximation is used, i.e., TSURF is a grid in z direction. The co-ordinate surfaces are given by z=TSURF(L)
(Default: NLTRZ = TRUE)
NLTRA = TRUE

toroidal approximation is used. The coordinate line is a polygonal approximation of a circle, or a segment thereof.
In case NLTOR = TRUE, the 3rd grid TSURF is a grid of toroidal angles. The toroidal segment (or the full torus) is approximated by NT3RD-1 straight cylinders.
In case NLTOR = FALSE, there are NTTRA toroidal periodicity boundaries, such that the toroidal segment is approximated, again, by NTTRAM=NTTRA-1 straight cylinders, but without toroidal resolution in the results. (Note: If NLTOR, NTTRA=NT3RD, internally.)
The torus axis of the entire mesh can be shifted in radial direction by adding a radial offset ROA (see below) to the x coordinates of the poloidal mesh (RSURF,PSURF).
The radial shift of the poloidal mesh RMTOR of the approximated torus is computed from ROA such that the volume inside radial surface NR1ST is exactly equal to the volume of an exact torus with poloidal cross section defined by the shifted radial grid (first standard mesh: RSURF).
Due to the approximations made by defining a torus by NTTRAM = NTTRA - 1 straight cylinders, this condition is fulfilled only approximately for the other radial surfaces. RMTOR converges to ROA with increasing NTTRA. NTTRA ≅ 30 is already a very good approximation.
(Default: NLTRA = FALSE).
NLTRT = TRUE

torus co-ordinates R,PHI,THETA. Presently being developed for NLSLB,NLPLG and NLTRI options. Not ready for use.
NT3RD
Number of grid-points in z- or toroidal direction
(default: NT3RD = 1, i.e. no grid is defined)
NTTRA
only needed in case NLTRA and .NOT.NLTOR. See above.
ZIA,ZGA,ZAA,ZZA,ROA

The 3rd grid ZSURF is defined in the same way as the x grid, using the parameters ZIA,ZGA,.... (cm) instead of RIA,RGA,....
In case of NLTRZ = TRUE , a z-grid is defined. ROA is irrelevant.
In case of NLTRA = TRUE , ZIA and ZAA are toroidal angles (in degrees). A grid of toroidal angles is defined. For example use ZIA=0.0 and ZAA=360.0 for a full torus. In case NLTOR this grid also defines the toroidal resolution. If .NOT.NLTOR, then periodicity at the endpoints ZIA and ZAA is automatically enforced.
(Defaults: NLTRA = FALSE, ZIA = 0 , ZAA = 1)
Data for mesh multiplication
NLMLT
the complete Standard mesh data are copied NBMLT times
NBMLT
Number of identical copies of the standard mesh. Transition from one such mesh (called "block" in EIRENE) into another one has to be defined by transparent additional surfaces (see block 3B)
VOLCOR
The volumes of all cells of the standard mesh as computed by EIRENE (in subroutine. VOLUME) are multiplied by one common factor VOLCOR for each "block".
Data for additional zones outside the Standard mesh
There are cells in the computational volume, which are defined through "additional surfaces" as cell boundaries. E.g. a standard mesh (is there is one) is augmented by "additional cells" in this case. If NLADD = FALSE , the complete block 2D may be omitted and the defaults are used.
Number of additional zones. The cell indexing along test flights in these zones has to be specified explicitly by making use of the ILSWCH, ILCELL parameters (block 3B)
Volume (cm−3) of each additional zone as seen by the test-particles.

2.2.1  Mesh Parameters

As stated above, the computational volume in an EIRENE run is subdivided into NSBOX cells,
NSBOX = (NR1ST · NP2ND · NT3RD) · NBMLT + NRADD .
The can mesh consist two different parts.
The first is the so called "standard mesh" part, defined by the grids RSURF, PSURF and TSURF. There may be NBMLT copies of that mesh, separated by additional surfaces (see below, section 2.3.2 (default: NBMLT=1).
The second part is the so called "additional cell region". Up to NRADD additional cells, of arbitrary shape, are defined by their boundaries: the "additional surfaces" defined in section 2.3.2.

NR1ST is the number of surfaces in the 1st (radial -or x-) grid. There are NR1STM = NR1ST - 1 cells.
The cell index of a particular cell is called NRCELL.
1 ≤ NRCELL ≤ NR1STM
NP2ND is the number of surfaces in the 2ND (poloidal -or y-) grid. There are NP2NDM = NP2ND - 1 cells.
The cell index of a particular cell is called NPCELL.
1 ≤ NPCELL ≤ NP2NDM
NT3RD is the number of surfaces in the 2RD (toroidal -or z-) grid. There are NT3RDM = NT3RD - 1 cells.
The cell index of a particular cell is called NTCELL.
1 ≤ NTCELL ≤ NT3RDM
For the cells of the standard mesh, NACELL = 0
For the cells in the "additional cell region" :
and
NRCELL = 0, NPCELL = 1, NTCELL = 1

The last radial (or x-) cell (number NR1ST), for any NPCELL, NTCELL inside the standard mesh, is not a real geometrical cell. It is used to store the radial (or x-) average of the volume averaged tallies.
Likewise, the last poloidal (or y-) cell (number NP2ND) is not a real geometrical cell, but used to store the poloidal (or y-) average of the volume averaged tallies (for any NRCELL, NTCELL inside the standard mesh)
NTCELL = NT3RD: analogically.

A particular cell can be specified in one of two possible ways:
Either by giving the 5 cell numbers:
NRCELL, NPCELL, NTCELL, NBLOCK, NACELL
or, alternatively, by giving the position in the 1-dimensional cell array
NCELL

The relation between these two options is:
NCELL = NRCELL + ((NPECLL-1)*NR1P2) + (NTCELL-1))*NP2T3 + NBLCKA
with
NBLCKA = (NBLOCK -1) * NSTRDT +NRADD

2.3  The Input Block for Surfaces

Grid surfaces may be assigned special properties (e.g.: reflecting, absorbing, periodicity, modified cell index switching, etc...). Their definition is described in subsection 2.3.1.
In addition to these grid surfaces also additional surfaces can be seen by the histories (vacuum boundaries, special surfaces for scoring fluxes, diagnostic surfaces, ....). Such surfaces can be general linear or second order surfaces in 3D space. Their definition and properties are described in subsection 2.3.2.

2.3.1  The Input Block for "Non-default Standard Surfaces"

General remarks

Some of the "standard grid surfaces" RSURF, PSURF and TSURF may be defined as reflecting or absorbing (rather than the default transparent character of co-ordinate surfaces), or as transparent but with other switching features (cell indexing, transition into additional cell region etc.) than the default settings. E.g., surface averaged tallies (see block 2.11, end of section 3.2 and Table 5.3) are only evaluated on either "Additional Surfaces" or on these "Non default Standard Surfaces", whereas they are not computed for those standard surfaces which are not explicitly identified in this block.
The geometrical properties of these "Non default Standard Surfaces" are specified here. This is described below.
Particle-surface interaction models
Although the (local) data for particle-surface interaction models (PSI-models) for each specific surface can be read in this input block, their meaning is described in block 2.6 together with the globally valid input data for particle-surface interactions.

*** 3A. Data for "Non-default Standard Surfaces"
 NSTSI
DO 31 ISTSI=1,NSTSI


unless otherwise stated, N = NLIM+ISTSI is the index of the following data:
    TXTSFL
ISTS IDIMP INUMP(ISTSI,IDIMP) IRPTA1 IRPTE1 IRPTA2 IRPTE2 IRPTA3 IRPTE3
ILIIN ILSIDE ILSWCH ILEQUI ILTOR ILCOL ILFIT ILCELL ILBOX ILPLG


optional for non transparent surfaces (ILIIN > 0):
The next 3 (or 4) lines comprise the block for local particle-surface interaction data (in analytic terminology: the boundary condition at this surface element). If they are omitted, the default particle-surface model is activated for this particular surface element, see section 2.6.
    ILREF  ILSPT ISRS  ISRC
ZNML EWALL EWBIN TRANSP(1,N) TRANSP(2,N) FSHEAT
RECYCF  RECYCT  RECPRM  EXPPL  EXPEL  EXPIL
RECYCS  RECYCC  SPTPRM  (this line may be omitted,
then: default sputter model, see \ref{sec2.6})
31 CONTINUE


also, optional for non-transparent surfaces (ILIIN > 0),
read a surface interaction model identifier to link one of the surface "local reflection models" defined in block 2.6 below to this current surface. The presence of such a link is identified by the string ′SURFMOD ′, followed by a name ′modname′. Later, in input block 2.6, a "local reflection model" block with that name ′modname′ must be included. This option allows quick changes of particle-surface interaction parameters, affecting many surfaces at once, by changes in just one location of the input file.
   SURFMOD_modname



Meaning of the Input Variables for "non-default standard surfaces"

NSTSI
Total number of non-default standard surfaces that do not act as prescribed by the default transparent standard co-ordinate surface model.
TXTSFL
Text to characterize a surface ("name of the surface") on the printout file.
ISTS
irrelevant; labelling index
IDIMP
flag to identify mesh from which this particular surface is chosen.
= 1
surface from the x- (radial) standard mesh (RSURF) Note that for the unstructured grid options NLTRI and NLTET (i.e. for the 2D triangular grid option or for the general 3D grids of tetrahedra) all surfaces are referred to as 1st grid (x or radial) surfaces, by abuse of language.
= 2
surface from the y- (poloidal) standard mesh (PSURF)
= 3
surface from the z- (toroidal) standard mesh (TSURF)
INUMP
Number of the surface in mesh RSURF, PSURF or TSURF respectively
IRPTA,IRPTE

Only a subregion of the surface acts by the "non-default options" specified for this particular surface. This subregion is defined by these flags.
If JMP is a surface from the first mesh, then IRPTA2→IRPTE2 and IRPTA3→IRPTE3 are the surface index ranges of the 2nd and 3rd mesh, respectively, for which this surface acts as non-default surface. IRPTA1 and IRPTE1 are irrelevant.
If JMP is a surface from the 2nd mesh, then IRPTA1→IRPTE1 and IRPTA3→IRPTE3 are the surface index ranges of the 1st and 3rd mesh, respectively, for which this surface acts as non-default surface. IRPTA2 and IRPTE2 are irrelevant.
If JMP is a surface from the 3rd mesh, then IRPTA1→ IRPTE1 and IRPTA2→ IRPTE2 are the surface index ranges of the 1st and 2rd mesh, respectively, for which this surface acts as non-default surface. IRPTA3 and IRPTE3 are irrelevant.
The third card ILIIN, ..., ILPLG is identical to the corresponding card for "additional surfaces", see block 2.3.2 below, "Input Data for Additional Surfaces". One exception is the flag ILSIDE, which controls the sign of the surface normal vector (hence the orientation of a surface). In case of unstructured "standard grids" NLTRI or NLTET (triangles in 2D and tetrahedrons in 3D) there is no well defined default surface orientation and the flag ILSIDE is irrelevant in such cases.

The input data in the block for the local reflection model are described below, see 2.6 : "Input Data for Particle-Surface Interaction Models".

2.3.2  Input Data for "Additional Surfaces"

General remarks

Internally each additional surface is defined by an algebraic equation and some algebraic inequalities specifying the boundary of that surface, i.e., that part of the surface which is seen by the test particles. By changing the sign of all coefficients in the algebraic equation the orientation of the surface normal vectors are reversed.
The intersection with the nearest surface along a trajectory is found by checking surfaces in the sequence of their input (in subroutine TIMEA of the geometry block GEO3D). This must be taken into consideration, if there are two parts of one surface specified by different boundaries and with non-empty intersection. In this case always the surface later in the input sequence is seen by the test particles.
One can define two distinct plane surfaces as one surface of second order, provided that this is compatible with the options available for surface boundaries.
We will refer to positive and negative directions of flights for particles intersecting a surface. By "positive" it is meant that the angle between the trajectory of the particle and the surface normal at the point of intersection is less than 90 degrees, "negative" has the opposite meaning.
speeding up geometry calculations:
Checking surfaces along a test flight can be abandoned depending upon the current position of a test particle.
For particles in cell no. NCELL all surfaces MSURF with IGJUM3(NCELL,MSURF)=1 are not checked. For particles on surface MS all other surfaces MSURF with IGJUM1(MS,MSURF)=1 are not checked. These matrices IGJUM1 and IGJUM3 can be set either in this block, see CH1 cards and CH3 cards described below, or in some user-specified problem specific routines (e.g. GEOUSR).
ignorable coordinates:
Surfaces with ignorable coordinates are specified either by setting the corresponding coefficients in the surface equation to zero, or, in case of 2-point, 3-point, 4-point of 5-point options, by setting the corresponding coordinates of the points to -1.D20 and +1.D20, respectively.
coordinate systems, transformations
All surfaces are specified in cartesian coordiantes.
In case of NLTRA (toroidal geometry approximation), surfaces must be defined in the local coordinate system of toroidal cells centered at the "magnetic axis", (xa,ya) = (RMTOR,0.), i.e., excluding the major radius of the torus.
To facilitate input of the geometrical data, each single (or a set of) surface(s) can be transformed by a Cartesian mapping after specification of the surface coefficients and boundaries in some "convenient" coordinate frame, by including "TRANSFORM cards" at the end of an input block for a particular surface. One such "TRANSFORM-deck" can act on an entire range of surfaces, but only on those which have been read previous to the "TRANSFORM-deck". Hence, in order to transform the entire set of additional surfaces by one deck, this deck must come after the last surface (no. NLIMI, see below).
Particle-surface interaction models
Although the (local) data for particle surface interaction models for each specific additional surface can be read in this input block, their meaning is described in block 2.6 together with the globally valid input data for particle surface interactions.

The Input Block

*** 3B. Data for Additional Surfaces
 NLIMI
CH0-cards" (may be omitted)
(format: 3A,69A)


arbitrary number of strings ±n/m, separated by blanks. n and m must be integer variables 1 ≤ n,mNLIMI.
 DO  ILIMI=1,NLIMI
TXTSFL
CH1-cards" (may be omitted)
(format: 3A,69A)


arbitrary number of strings ±n/m, separated by blanks. n and m must be integer variables 1 ≤ n,mNLIMI.
   CH2-cards" (may be omitted)
(format: 3A,69A)


arbitrary number of strings ± n/m, separated by blanks. n and m must be integer variables 1 ≤ n,mNLIMI.
   RLBND  RLARE  RLWMN  RLWMX
ILIIN ILSIDE ILSWCH ILEQUI ILTOR ILCOL ILFIT ILCELL ILBOX ILPLG
IF (RLBND.LT.0.) THEN
A0LM  A1LM  A2LM  A3LM  A4LM  A5LM
A6LM  A7LM  A8LM  A9LM


let RLBND = -KL, then first read L cards:
     ALIMS  XLIMS  YLIMS  ZLIMS


and then read K blocks next:
     ALIMS0  XLIMS1  YLIMS1  ZLIMS1  XLIMS2  YLIMS2
ZLIMS2  XLIMS3  YLIMS3  ZLIMS3
ENDIF
IF (RLBND.EQ.0.) THEN
A0LM  A1LM  A2LM  A3LM  A4LM  A5LM
A6LM  A7LM  A8LM  A9LM
ENDIF
IF (RLBND.GT.0..AND.RLBND.LT.2.) THEN
A0LM  A1LM  A2LM  A3LM  A4LM  A5LM
A6LM  A7LM  A8LM  A9LM
XLIMS1  YLIMS1  ZLIMS1  XLIMS2  YLIMS2  ZLIMS2
ENDIF
IF (RLBND.GE.2.) THEN
P1(1,..)  P1(2,..)  P1(3,..)  P2(1,..)  P2(2,..)  P2(3,..)
IF (K .GT. 2) THEN
P3(1,..)  P3(2,..)  P3(3,..)  P4(1,..)  P4(2,..)  P4(3,..)
IF (K .GT. 4) THEN
P5(1,..)  P5(2,..)  P5(3,..)
ENDIF
ENDIF
ENDIF


optional for non transparent surfaces (ILIIN > 0):
The next 3 (or 4) lines comprise the block for local particle-surface interaction data (in analytic terminology: the boundary condition at this surface element). If they are omitted, the default particle-surface model is activated for this particular surface element, see section 2.6.
    ILREF  ILSPT ISRS  ISRC
ZNML EWALL EWBIN TRANSP(1,N) TRANSP(2,N) FSHEAT
RECYCF  RECYCT  RECPRM  EXPPL  EXPEL  EXPIL
RECYCS  RECYCC  SPTPRM  (this line may be omitted,
then: default sputter model,
see \ref{sec2.6})


also optional for non transparent surfaces (ILIIN > 0):
read a surface interaction model identifier to link one of the surface "local reflection models" defined in block 2.6 below to this current surface. The presence of such a link is identified by the string ′SURFMOD ′, followed by a name ′modname′. Later, in input block 2.6, a "local reflection model" block with that name ′modname′ must be included. This option allows quick changes of particle-surface interaction parameters, affecting many surfaces at once, by changes in just one location of the input file.
   SURFMOD_modname


optional: read one or several blocks of five cards each for orthogonal mapping. The presence of each such block is identified by the first card of that block containing the string ′TRANSFORM′ followed by the other four cards:
   TRANSFORM
ITINI,ITEND
XLCOR,YLCOR,ZLCOR
XLREF,YLREF,ZLREF
XLROT,YLROT,ZLROT,ALROT


 ENDDO   (end of ILIMI loop)



Meaning of the Input Variables for additional surfaces

CH0 n1/m1 n2/m2 ...
surfaces from the range n1 to m1, n2 to m2, ..., are ignored by EIRENE. Specifying a surface in such a CH0-card is identical to taking it out from the input file. It may be more convenient in some cases, however, to use the CH0 option, because the labelling index of the remaining valid surfaces is not altered then. No input is read for these surfaces, and the input segment for the next valid surface (identified by the string '*text') is read next.
NLIMI
number of surfaces in the input block
TXTSUR
Text to identify a surface ("name of the surface") on the printout file
CH1(ILIMI) n1/m1 n2/m2 ...
surfaces from the range n1 to m1, n2 to m2, ..., are considered invisible for a particle located on this current surface ILIMI. Intersection of trajectories starting from surface ILIMI with those "invisible" surfaces is not checked.
CH2(ILIMI) n1/m1 n2/m2 ...
Only for second order surfaces ILIMI, n1-m1, ... The first of the two possible intersections is ignored for particles located on surface no. ILIMI.
General Data for Surfaces
RLBND
flag for different options to define the boundary of the surface
RLARE
Area (in cm2) of the surface element which is seen by the test particles. (Default: 666.0) (needed only for scaling of non default surface averaged tallies) If RLARE is not specified here, (i.e., if a value less than or equal to zero is read) then EIRENE tries to evaluate this area itself. For some surfaces this is still not possible automatically.
RLWMN
lower weight limit for space weight window for particles crossing the surface in positive direction. (not in use)
RLWMX
upper weight limit for space weight window for particles crossing the surface in positive direction. (not in use)
ILIIN
defines the type of surface
> 0
non-transparent surface
= 1
reflecting, partly or purely absorbing surface.
local reflection model has to be specified unless default model is to be used; all surface tallies (see Table 5.3) are updated and a switch can be operated.
= 2
purely absorbing surface ;. surface tallies for incident fluxes (i.e., POT... and EOT... tallies in Table 5.3) are updated and the particle history is stopped then.
= 3
mirror for incident test particles. I.e., specular reflection for neutral test particles, and for charged test particles the sign of the velocity component parallel to the B-field is reversed.
=m4
periodicity surface, with regard to x, y, or z coordinate, depending upon whether this surface is a standard x, y, or z grid surface, respectively. Move particle to x / radial grid surface no. m, m integer (or to y / poloidal or to z / toroidal surface no. m, respectively) and continue track from there with otherwise identical particle parameters.
This option is currently implemented only for Cartesian grids (NLSLB and NLTRZ) and for non-default standard grid surfaces only.
≤ 0
transparent surface (for example: hole in one of the other "additional surfaces"). Particle and energy fluxes onto and from these surfaces do not contribute to global balances.
= 0
Particle history is not interrupted
No surface tallies are updated, no switches can be operated. Fastest option.
= -1
Particle history will be stopped and restarted. A switch can be operated. I.e., this surface is used only for switching (see below: ILSWCH) or re-initializing the particle's track at the point of intersection. No surface tallies are updated.
= -2
if a particle is crossing the surface in the positive direction, (one sided-) surface tallies are updated, e.g., by default: partial particle and energy currents J+ (Amp) and K+ (Watt). These are stored in the POT... and EOT... tallies of Table 5.3. If the particle crosses the surface in the direction opposite to the surface normal, then negative partial particle and energy currents J (Amp) and K (Watt) are updated. These are stored in the PRF... and ERF... tallies of Table 5.3.
= -3
Net currents (e.g. J+J), are evaluated, and stored on the POT... and EOT... tallies (see Table 5.3). The PRF... and ERF... tallies are empty for these surfaces.
≤ -4
Not in use. Currently: same as ILIIN=-2 option.
ILSIDE

= 0
both sides of the surface act as described by ILIIN option (default).
= 1
particles incident on the surface in the negative direction will be absorbed (i.e., ILIIN = 2 option from that side).
= 2
particles incident on the surface in the negative direction will be killed and the message
or
"ERROR IN STDCOL"
will be printed. The contribution of these particles to the particle- and energy flux balances will be called PTRASH and ETRASH respectively. This option should be used for geometry testing whenever the user expects particles incident only from one side.
= 3
particles incident on the surface in the negative direction will not see the surface, i.e., this surface acts like a (semi) transparent surface (ILIIN = 0 option) from that side.
= -1
as 1, but with the opposite direction of the surface normal
= -2
as 2, but with the opposite direction of the surface normal
= -3
as 3, but with the opposite direction of the surface normal
ILSWCH
= IJKLMN, i.e. six digits I, J, K, L, M and N
= 0
no switch is operated
N
EIRENE flag ITIME
N = 1
The calculation of the step sizes in the standard mesh is abandoned for a particle which crosses the surface in the positive direction, and is reactivated, if the particle strikes in the negative direction
N = 2
as 1, but with the direction of the surface normal reversed for this option.
M
EIRENE flag IFPATH
M = 1
Abandon the calculation of the collision rates (entry into the vacuum) for a particle which is striking the surface in the direction of the surface normal. For particles incident from the other direction, evaluation of collision rates is reactivated.
M = 2
as 1, but with the direction of the surface normal reversed.
L
EIRENE flag IUPDTE
L = 1
Abandon the updating of volume-averaged tallies for a particle which is striking the surface in the direction of the surface normal. For particles incident from the other direction, updating of volume averaged tallies is reactivated.
L = 2
as 1, but with the direction of the surface normal reversed for this option.
I,J,K flags for switching cell numbers at transition into a different mesh cell.
K
for particles in an additional cell, i.e., not in one of the "standard mesh" blocks:
K = 1
Increase the actual additional cell number NACELL for a particle striking the surface in the direction of the surface normal by ILACLL. Decrease NACELL by ILACLL if the particle is striking in the negative direction. Specification of ILACLL is via the input variable ILCELL, see below.
K = 2
as K = 1, but with the direction of the surface normal reversed for this option.
for particles inside the "standard mesh ", i.e., not in the "additional cell region"
K = 1
Increase the standard mesh block number NBLOCK for a particle striking the surface in the direction of the surface normal by ILBLCK. Decrease NBLOCK by ILBLCK if the particle is striking in the negative direction. Specification of ILBLCK is via the input variable ILCELL, see below.
K = 2
as K = 1, but with the direction of the surface normal reversed.
J
for particles at the boundary between "additional" and "standard" mesh regions.
J = 1
entrance into standard mesh, block no.
NBLOCK = ILBLCK
or exit from standard mesh into additional cell
NACELL = ILACLL.
Specification of ILACLL and ILBLCK is via the input variable ILCELL, see below. If ILACLL = 0, then no switch to additional cell is operated. (E.g.: for surfaces which are reflecting from this side).
J = 2
as J = 1. The direction of the surface normal does not matter here.
I
similar to J-flag, i.e., for transitions between standard and additional meshes, but different cell number switching.
I = 1
Entrance into standard mesh, block no.
NBLOCK = NACELL+ILBLCK,
if the particle is striking in the positive direction, or
NBLOCK = NACELL-ILBLCK,
if the particle is striking in the negative direction. Exit from standard mesh into additional cell
NACELL = NBLOCK+ILACLL,
for a particle striking the surface in the positive direction, or
NACELL = NBLOCK-ILACLL,
if the particle is striking in the negative direction.
Specification of ILACLL and ILBLCK is via the input variable ILCELL, see below.
I = 2
as I = 1, but with the direction of the surface normal reversed.
If a test particle history starts from a surface (NLSRF option), then ILSWCH acts as if this particle had struck the surface prior to the birth process in the positive direction. This default setting is only available for ILSIDE ≠ 0 and can (must) be overruled by the SORIFL flag (see section 2.7), e.g. if a surface source needs to be defined on a surface with ILSIDE = 0.
ILEQUI
The algebraic equations for the surfaces J and IABS(ILEQUI(J)) will be described by exactly the same coefficients (up to a common sign, if ILEQUI(J) .lt. 0). For example a triangle can be specified by the three corners and another part of the same plane surface can be specified directly by its algebraic coefficients. To avoid round-off errors one should use the ILEQUI option in such cases, in particular if surface J is a transparent "hole" in surface ILEQUI(J), or vice versa.
Default: ILEQUI = 0
ILTOR
For NLTRA option only (see block 2c):
if ILTOR > 0 :
the surface is defined with respect to the local coordinate system of the toroidal segment with "cell-number" ILTOR, hence: 1 ≤ ILTOR ≤ NTTRAM.
if ILTOR=0 :
the surface is defined with respect to any local coordinate system. I.e., the surface equations are taken to be the same in each local system.
If the surface equations are z-independent, then this surface is toroidally symmetric (within the NLTRA-approximation).
Otherwise the surface has NTTRAM-fold periodicity.
This flag is irrelevant for NLTRZ (i.e., if cylindrical coordinates are used) or for NLTRT.
Default: ILTOR = 0
ILCOL
Flag for the color that is used for plotting this surface on 2d or 3d geometry plots. If ILCOL ≤ 0 than -ILCOL is used and the surface area is filled in by that color on the 3d geometry plots.
Default: ILCOL=1
ILFIT
= MN
This option is relevant only for surfaces with one ignorable coordinate, i.e. it only works for the 2. ≤ RLBND ≤ 3. surface boundary options.
It is a tool to facilitate a neat fitting of surfaces, in particular for connecting curved and plane surfaces (i.e. to avoid particle leakage due to numerical round-off errors in the algebraic surface coefficients.
M and N (3 digits each, M may be omitted if not needed) are the numbers of the surfaces (which must have the same ignorable coordinate) the boundaries of which should match to those of the actual surface J. The boundaries of these neighboring surfaces M and N must be specified by the RLBND = 1 or RLBND = 1.5 option.
The fitting is achieved by a small automatic internal modification of the surface data P1 and/or P2 of surface number J in subroutine SETFIT. Printout of the modifications made there is activated with the TRCSUR flag (input block 11).
ILCELL
Parameter ILBLCK and ILACLL for the ILSWCH flags described above. Let ILCELL = NM, with N and M being integers with 3 digits each. Then N = ILBLCK and M = ILACLL.
ILBOX
ILPLG EIRENE can write out information for a finite element mesh generator to produce a grid of triangles for a muliply connected 2D domain with cracks and holes. The various (inner and outer) boundaries are given as polygonal lines, which are composed of selected standard grid surface segments (NLPLG option) and/or additional surfaces (2 ≤ RLBND < 3 option).
This flag identifies closed polygonal lines composed of additional surfaces given by the 2-point option and/or of standard surfaces in the x-y-plane. For example if ILPLG(I)=NN, for surfaces I = I1, I2, ...IN, (NN a positive integer) then these N surfaces form a closed polygonal line in the x-y-plane. The region inside this closed line is part of the computational domain. By a negative integer value of NN a closed polygonal region can be excluded from the computational domain, i.e., a hole in the domain is specified by these surfaces. EIRENE writes an output file appropriate for a finite element mesh generator (available from FZ-Juelich) to produce a triangular discretization of the resulting (possibly multiply connected) domain. This option can be used to discretize arbitrarily complex 2D domains with internal and external boundaries given by the additional or non-default standard surfaces.
These finite element grids can be combined with the regular grids by using the problem specific geometry routines (see section 3) or the code interfacing routines INFCOP (see section 4).

Surface coefficients and boundaries of surfaces

The surface equation:
 A0LM
 +
 A1LM ·x + A2LM ·y + A3LM ·z +
 +
 A4LM ·x2 + A5LM ·y2 + A6LM ·z2 +
 +
 A7LM ·xy + A8LM ·xz + A9LM ·yz = 0
(3.1)
The positive surface normal (nx , ny , nz) depends upon the point of impact (x,y,z) and is defined by the vector
 nx
 =
 A1LM + 2 ·A4LM ·x +A7LM ·y + A8LM ·z
 ny
 =
 A2LM +A7LM ·x + 2 ·A5LM ·y +A9LM ·z
 nz
 =
 A3LM +A8LM ·x + A9LM ·y +2 ·A6LM ·z
For a plane surface the following reduction is valid:
(nx , ny , nz) = (A1LM,A2LM,A3LM)

The boundary of the valid part of the surface may be described by 4 different options, depending upon the value of the flag RLBND.
RLBND = 0

No boundary inequalities specified, i.e. the whole surface is seen by the test particles.
0 < RLBND < 2

= 1
Only that part of the surface, which lies inside the right parallelepiped defined by the two vectors (XLIMS1, YLIMS1, ZLIMS1) and (XLIMS2, YLIMS2, ZLIMS2), is seen by the particles. I.e. the three inequalities
XLIMS1 ≤ xXLIMS2
YLIMS1 ≤ yYLIMS2
ZLIMS1 ≤ zZLIMS2
are checked at the point of intersection (x,y,z).
= 1.5
Complement to RLBND = 1.
Only the surface element outside the parallelepiped is seen by the particles.
RLBND ≥ 2

In this case the surface will be defined by the input of the coordinates of at least 2 and at highest 5 points on a plane surface. If there are only 2 points, the surface is parallel to one axis. If there are 3 or more points, then the boundary of this plane surface is a closed polygon (P1 ,..., Pn , P1). Therefore, the correct order of points at input is relevant. The orientation of the positive surface normal vector is defined by the first 3 points, and it is given by the vector product (P3 − P1) ×(P3 − P2). Thus, the orientation can be reversed e.g. by interchanging P2 and P3.
2.1
plane surface parallel to z axis. The surface equation of this plane reads ax+by+c=0 with the coefficients a,b and c such that the points P1, P2 lie on this surface and the valid part of that surface ranges from P1 to P2 in the xy-plane. The z-coordinates of these two points define the boundaries in z direction
= 2.2
Complement to RLBND = 2.1
= 2.4
as RLBND=2.1 option, but with z and y exchanged. I.e., now the y coordinates of the points P1, P2 are the boundaries of the surface ax+bz+c=0 in y direction.
= 2.5
Complement to RLBND = 2.4
= 2.7
as RLBND=2.1 option, but with z and x exchanged. I.e., now the x coordinates of the points P1, P2 are the boundaries of the surface ay+bz+c+0 in x direction.
= 2.8
Complement to RLBND = 2.7
= 3
plane triangle defined by the corners P1 ,P2 , P3
P1=(P1(1),P1(2),P1(3))
P2=(P2(1),P2(2),P2(3))
P3=(P3(1),P3(2),P3(3))
= 3.5
complement to RLBND = 3; only The plane surface outside the triangle is seen by the test particles.
= 4
plane quadrangle; surface inside the polygon
(P1 , P2 , P4 , P3 , P1).
Here P1, P2, P3 are as in the RLBND=3 option, and P4 = (P4(1), P4(2), P4(3))
Thus this surface is the union of the triangles with vertices P1, P2, P3 and P2, P4, P3 respectively.
= 4.5
complement to RLBND = 4; only the part of the plane surface outside the quadrangle is seen by the test particles
= 5
plane quint-angle; surface inside the polygon (P1, P2, P4, P5, P3, P1) P1, P2, P3, P4 as RLBND=4, and P5 = (P5(1), P5(2), P5(3))
= 5.5
complement to RLBND = 5; only the part of the plane surface outside the quint-angle is seen by the test particles.
RLBND < 0

-KL
The surface is bounded by L linear inequalities and by K second order inequalities.
ALIMS + XLIMS * x + YLIMS * y + ZLIMS * z ≤ 0
(L inequalities)

 ALIMS0
 +
 XLIMS1 ·x + YLIMS1 ·y + ZLIMS1 ·z
 +
 XLIMS2 ·x2 + YLIMS2 ·y2 + ZLIMS2 ·z2
 +
 XLIMS3 ·xy + YLIMS3 ·xz + ZLIMS3 ·yz ≤ 0
(K inequalities)
The meanings of the input variables for the local reflection model are described below (see: 2.6: "Input Data for Surface Interaction Models").
The meanings of the input variables for a transformation block are:
ITINI, ITEND

The transformation defined by the next 3 cards is carried out for additional surfaces from number ITINI through ITEND The transformation is carried out as soon as this transformation deck is found. Hence, if such a transformation deck is found after "additional surface" no. IS, then one must guarantee: ITINI ≤ IS and ITEND ≤ IS.
XLCOR, YLCOR, ZLCOR

Translation:
The origin of the coordinate system is shifted by the vector
XLCOR,YLCOR,ZLCOR
XLREF, YLREF, ZLREF

Reflection:
to be written
No transformation if XLREF,YLREF,ZLREF = (0,0,0)
XLREF, YLREF, ZLREF

Rotation:
The vector XLROT,YLROT,ZLROT defines the axis of rotation.
ALROT is the angle of rotation (in degrees)
No transformation if ALROT=0 or axis of rotation= (0,0,0)

2.4  Input Data for Species Specification and Atomic Physics Module

General remarks

EIRENE can handle up to NATM "atom" species, NMOL "molecule" species, NION "test ion" species, NPHOT "photon" species (lines) and NREAC different atomic, molecular or photonic reactions between these "test particles" and the "bulk ions" or electrons. There may be up to NPLS "bulk ion" species, and one electron gas derived internally from the assumption of local charge neutrality. Amongst the heavy background particles (the "bulk ions") may be species with charge state zero, i.e., neutral particles. By abuse of language, we will refer to them as "bulk ions" as well, but meaning in this case: heavy background particles.
NATM, NMOL, NION, NPHOT, NPLS and NREAC are specified in the parameter block "PARMUSR", see: "Parameter Statements" (section 3.1 below, for versions EIRENE2000 and older), or are determined automatically from the input file.
Masses in EIRENE (RMASSA, RMASSM, RMASSI, RMASSP for atoms, molecules, test ions and bulk ions, respectively) are in units of the proton mass, which is taken to be 1.0073 in atomic mass units.
First EIRENE reads "non-default" atomic data (if any, i.e., if NREAC > 0) specified in so called "reaction cards", IR ≤ NREACI ≤ NREAC

IR  FILNAM  H123  REAC  CRC  MASSP  MASST  DP RMN  RMX

(format: I3,1X,A6,1X,A4,Axxx,A3,2I3,3E12.4)


from external A&M data files such as FILNAM = HYDHEL, = METHAN, or = AMJUEL, or other data files of appropriate format.
For collision reaction rates, electron cooling rates, etc., from the "ADAS database format AFD11", or any other properly tabulated data source, i.e. for FILNAM = ADAS, an additional card is read for each IR entry, see below, to identify the species and ion charge state within the ADF11 file.

IR  ADAS  H123  REAC  CRC  MASSP  MASST  DP

(format: I3,1X,A6,1X,A4,Axxx,A3,2I3,E12.4)
ELNAME, IZ

(format: 4X,A2,1X,I3)


For photonic emission and absorption data (line shapes), i.e. for FILNAM = PHOTON, the line is shortened to

IR  FILNAM  H123  REAC

(format: I3,1X,A6,1X,A4,Axxx)


for reading from the spectroscopic line-shape database PHOTON, but additional lines are needed to specify the line broadening options. See below.
There also exist an option for specifying constant (or very simple functional dependencies) cross sections, reaction rate coefficients, reaction rates, etc. directly via the input file, without resorting to an external database. This has proven useful for testing purposes, e.g. comparison with simple analytic or numerical solutions. In this case: FILNAM = CONST, further details: see below.
Later EIRENE will assign the relevant atomic molecular of photonic processes to each test particle (and also bulk particle) from either this set of NREACI processes, or from a small set of "default processes", which are hard wired in EIRENE.
In versions younger than 2002 each reaction in these databases can contain a string 'fit-flag = value', with IFTFLG=value. IFTFLG is used in EIRENE to identify the type of fitting expression to be evaluated with the fitting coefficients from the database.
By default: IFTFLG=0 for all data, i.e., single or double polynomial fits.
All atomic rate coefficients and cross sections can be scaled with a constant factor (FREAC, see below), for sensitivity studies.
Excluded from this scaling option are all photonic processes and those elastic collisions for which EIRENE uses an interaction potential (H.0-type data), because here the scaling would not have the expected effect on transport, but rather it would only modify the effective small angle scattering cut-off in the binary collisions.
Atomic and molecular processes are always of the following type:

 a ·A + b ·B → m ·M + n ·N + ∆E
Here A,B,M and N are labels for the type of pre- and post collision particles. a,b,m and n are the stoichiometric coefficients, and ∆E is the amount of internal energy transferred into (or at the expense of) the kinetic energy of the collision products M and N.
The following conventions are always in use:
At least one of A, B, M or N must stand for a "test particle", for otherwise this process is not relevant for the transport problem solved by EIRENE.
• particles A are always "bulk particles", i.e., from the list of "bulk ions" (input block 5) or electrons.
• particle B may be one "test particle". Then this process A + B→ ... must be included in the list of reaction decks for this test particle B. Hence: b = 1 then.
• particle B may be one "bulk particle". Then this process must be included in input block 7 as a volume-source for particles M and/or N,
• particles M and/or N can be either secondary "test particles" or secondary "bulk particles".
• non-linear processes, i.e., processes in which both A and B should, in principle, stand for test particles, can be included by specifying particle A as bulk particle in input block 5 and, simultaneously, as test particle in input block 4. The parameters for the "artificial bulk species" A are iterated in a sequence of EIRENE runs. In cases in which species M and/or N are also available as ("real" or "artificial") bulk species, there is the choice to specify secondaries from a reaction either as test-particles and to continue the history with those, or, alternatively, to specify them as bulk particles, and stop the trajectory. In this latter case one must add a volume source to launch these "secondaries" via a proper spatially distributed source in the next iteration. By using this option carefully, and combining it with the stratified source sampling technique, coupling to external codes/models can be made either more or less implicit.
Some care is needed here to avoid double-counting.
Next EIRENE expects one so called "species block" for each test particle species. Later, in input block 5, it will also expect one species block for each of the heavy background particles, i.e., for the "bulk ions".
A "species block" has the format:
ISPZ$TEXT$ NMASS$NCHAR$ NPRT$NCHRG$ ISRF$ISRT$ ID1$NRC$  NFOL$NGEN$  NHSTS$ID3$
(format: I2,1X,A8,1X,12(I2,1X))


In case ID1 ≠ 3, i.e., two or less post collision particle (not counting electrons)
 DO  K= 1, NRC$IREAC$  IBULK$ISCD1$  ISCD2$ISCDE$  IESTM$IBGK$
EELEC$EBULK$  ESCD1$EDUMMY FREAC$  FLDLM$ENDDO  New option in 2006: In case ID1 = 3, i.e., three post collision particles (not counting electrons)  DO K= 1, NRC$
IREAC$IBULK$  ISCD1$ISCD2$  ISCD3$ISCDE$  IESTM$IBGK$
EELEC$EBULK$  ESCD1$EDUMMY FREAC$  FLDLM$ENDDO  New option in 2006: In case ID1 = 4, i.e., four post collision particles (not counting electrons)  DO K= 1, NRC$
IREAC$IBULK$  ISCD1$ISCD2$  ISCD3$ISCD4$ ISCDE$IESTM$  IBGK$EELEC$  EBULK$ESCD1$  EDUMMY FREAC$FLDLM$
ENDDO


where $stands for either A (atoms), M (molecules), I (test ions), PH (photons) or P (bulk ions). I.e., a species block consists of one species specification card ISPZ$ ..., and NRC$"reaction decks" of two cards each. For some particle species (in particular for hydrogenic atoms, molecules and molecular ions, and for helium atoms) EIRENE has default A&M data, to which it resorts, if no reaction cards are in the input deck (i.e., if NRC$=0, see below) for a particular species. These default models are described at the end of this section. They are overruled by the reaction cards.
In order to de-activate all possible reactions (also the default reactions), e.g., to simulate the collision-less "Knudsen flow" for all or a particular test particle species, one must set NRC$=−1. Usually, a particular particle type and species in the species blocks is identified by an integer IPART$=LMN (3 digits). Here N fixes the type of the particle, M is the number of particles characterized by IPART$and L is the species index within the specified group (type) of particles. The following values of N are currently in use: • N=0: photons (versions 2004 and younger) • N=1: atoms • N=2: molecules • N=3: test ions • N=4: bulk ions • N=5: electrons The Input block *** 4. Data for Species Specification and Atomic Physics Module Note: in EIRENE2003 or later the format of input flag REAC described below was generalized from A9 to Axxx, with xxx ≤ 50, to accommodate the longer reaction information for photonic processes into the same format. The old format A9 is of course automatically still recognized as special case. NREACI DO  for "real particles", databases: FILNAM = HYDHEL, AMJUEL, METHAN, H2VIBR  IR FILNAM H123 REAC CRC MASSP MASST DP RMN RMX (format: I3,1X,A6,1X,A4,Axxx,A3,2I3,3E12.4) IF (RMN.GT.0.) READ IFEXMN, FPARM(J),J=1,3 IF (RMX.GT.0.) READ IFEXMX, FPARM(J),J=4,6  for "real particles", database: FILNAM = ADAS (or any equivalent ADF11 format)  IR ADAS H123 REAC CRC MASSP MASST DP (format: I3,1X,A6,1X,A4,Axxx,A3,2I3,3E12.4) ELNAME, IZ (format: 4X,A2,1X,I3)  i.e.: one additional input card is read to identify the chemical element and charge state. Extrapolation flags RMN and RMX are not in use in this case and may be omitted. Default (hard wired) extrapolation schemes are used if the plasma density and temperature is outside the tabulated range. for "real particles", database: FILNAM = CONST. Fit parameters are directly read from input file, up to nine coefficients. In EIRENE2005 or later: Optional : Further parameter FTFLAG, length: not more than 9 characters.  IR CONST H123 FTFLAG(optional) CRC MASSP MASST DP (format: I3,1X,A6,1X,A4,(A9),xxxX,A3,2I3,3E12.4)  i.e.: Reading of FTFLAG is optional. Default for FTFLAG: =0. Parameter REAC does not exist. Extrapolation flags RMN and RMX are not in use in this case and may be omitted. for "photons", database FILNAM = PHOTON  IR PHOTON H123 REAC (format: I3,1X,A6,1X,A4,Axxx) IPRFTYPE IPLSC3 IMESS IFREMD NRJPRT DO IFREMD II KENN IK6 (format: I6,1X,A2,3X,I6) ENDDO  I.e.: after the first line there is a special format for this set of input cards in case of photonic processes, read from the spectral database PHOTON. In case IFREMD > 0 the additional cards specify options for "foreign pressure line broadening". Then the subroutine SLREAC is called, which picks the atomic data (fit coefficients) for reaction "IR" from the file FILNAM. It then stores this information on the arrays in Module COMAMF with label IR.  * 4A. atoms species cards NATMI DO 44 IATM= 1, NATMI  * read NATMI species blocks with$ = A
 44 CONTINUE
* 4B. molecules species cards
NMOLI
DO 46  IMOL= 1, NMOLI


* read NMOLI species blocks with $= M  46 CONTINUE * 4C. test ions species cards NIONI DO 48 IION= 1, NIONI  * read NIONI species blocks with$ = I
 48 CONTINUE

and, for versions younger than 2003:

* 4D. test ions species cards
NPHOTI
DO 49  IPHOT= 1, NPHOTI


* read NPHOTI species blocks with $= Ph  49 CONTINUE  Meaning of the input variables NREACI Total number of different reactions to be read. The next block has different meanings for "real particles" and "photons". Cross section and rate coefficients are specified for particles, but emission and absorbtion line shapes are specified for photons. "real particles" An interaction potential V(r), cross section σ plus a reaction rate coefficient 〈σv〉 for one process counts as one reaction, but higher order rate coefficients such as energy or momentum weighted rate coefficients for the same process count as new reaction and must be labelled by a different index IR (see below). Storage is provided for up to NREAC different additional reactions (see "Parameter Statements", section 3.1), i.e., one must guarantee NREACI.LE.NREAC IR Labeling index of the reaction that is read in. The condition IR.LE.NREACI must be fulfilled (otherwise: error exit). FILNAM Name of the data-set that contains the required reaction. HYDHEL Hydrogen and Helium data [23] METHAN Hydrocarbons data [24] H2VIBR Data for H2 molecules and their isotopomeres, including effects of vibrational and electronic excitation, as well as rates needed for radiation trapping calculations. AMJUEL supplement to HYDHEL and METHAN for neutral gas transport Monte Carlo codes, e.g. multi-step reaction rates, etc. CONST reaction IR is a collision process with explicitly specified fit coefficients for cross sections or reaction rate coefficients, (depending upon input flag "H123"), e.g., constant, power law, etc. These fitting coefficients are directly read from the formatted input file rather than from an external file. ADAS Collisional radiative rate coefficients from ADAS ADF11 files, which must be located in a directory, the path to which is specified in input block 1, section 2.1, by one of the "CFILE" cards described there: CFILE ADAS path/adf11/ PHOTON H123 Identification flag for the type of data: interaction potential parameters, cross section, rate coefficient, momentum loss rate coefficient, energy loss rate coefficient, etc. case 0: FILNAM = CONST currently available only for: H123 = H.1, = H.2, = H.5, or = H.8, i.e., only for single parameter fits. in case IFTFLG = 10, 110, 210,... , "single parameter", constant collision parameter. One more line F(0), REAL is read. For all other values of IFTFLG, including the default IFTFLG=0, two more input lines with 9 fitting coefficients F(I), I=0,5 (REAL) F(I), I=6,8 (REAL) are read, in subroutine SLREAC from the formatted input file. The natural logarithm of cross-section: ln(σ) (H123 = H.1) or of a rate coefficient: ln(R) (H123 = H.2, = H.5, or = H.8) is computed as: ln(σ) = ∑I=08 F(I) ·ln(Elab)I and, likewise, a rate coefficient R is evaluated as: ln(R) = ∑I=08 F(I) ·ln(T)I Here Elab is the relative energy of impact, but with the mass of the charged particle (i.e., for a target of neutral particles at rest, see below), and T is the electron or ion temperature, depending on the type of impacting plasma particle. This option permits specification of constant cross sections or rate coefficients (only the first of the nine parameters nonzero), or, e.g., specific power law energy dependencies or temperature dependencies (second of the nine parameters). case 1: FILNAM = HYDHEL, METHAN, AMJUEL, H2VIBR H.1 single parameter fit for cross section σ(cm2) as function of energy ELAB (eV), with ELAB = mLab / 2 ·vrel2 (For the definition of mLab = MASSP : see below. H.2 single parameter fit for rate coefficient 〈σv 〉(cm3/s) as function of target temperature (eV), assuming "zero velocity" projectile H.3 double parameter fit of rate coefficient 〈σv〉(cm3/s) as function of projectile energy (eV) and target temperature (eV) H.4 double parameter fit of rate coefficient 〈σv 〉(cm3/s) as function of target density (cm−3) and target temperature (eV) H.5 single parameter fit of target particle momentum weighted rate coefficient 〈 σ v · mp · |vp| 〉 (cm3/s ·AMU ·cm/s) as function of target temperature (eV) H.6 double parameter fit of target particle velocity weighted rate coefficient 〈σv ·mp ·| vp | 〉(cm3/s ·AMU ·cm/s) as function of projectile energy (eV) and target temperature (eV) H.7 double parameter fit of target particle velocity weighted rate coefficient 〈σv ·mp ·| vp | 〉(cm3/s ·AMU·cm/s) as function of target density (cm3) and target temperature (eV) H.8 single parameter fit of target particle energy weighted rate coefficient 〈σv ·mp /2 ·vp2 〉(cm3/s ·eV) as function of target temperature (eV) H.9 double parameter fit of target particle energy weighted rate coefficient 〈σv ·mp /2 ·vp2〉(cm3/s ·eV) as function of projectile energy (eV) and target temperature (eV) H.10 double parameter fit of target particle energy weighted rate coefficient 〈σ v · mp/2 · vp2〉 (cm3/s ·eV) as function of target density (cm3) and target temperature (eV) H.11 single parameter fit for any other data, e.g. to be used in special user supplied programs. H.12 double parameter fit for any other data, e.g. to be used in special user supplied programs, (i.e. not understood by EIRENE) case 2: FILNAM = ADAS H.4 double parameter table of rate coefficient 〈σv 〉(cm3/s) as function of target density (cm−3) and target temperature (eV) read from ADAS adf11 files. Filenames starting with scd... (ionization rate coefficients contain ionization rate coefficients to charge state Z, Z=1,... Zmax) and filenames starting with acd...... contain recombination rate coefficients, from charge state Z, Z=1, Zmax. The particular value of Z is specified in the next input card. case 3: FILNAM = PHOTON to be written REAC Name or number of the reaction as used in the file FILNAM (e.g.: REAC = 2.15.2 in file METHAN for the reaction e + CH4+CH3+ + H + e ). If FILNAM = CONST, then REAC can be used to input the flag IFTFLG. Whether or not REAC is to be interpreted as IFTFLG is identified by presence of the string FT (e.g.: 'REAC' = 'FT 110' would lead to IFTFLG = 110 for this reaction). CRC Identification for the type of the reaction process. Depending upon the value of this flag various assumptions are automatically made concerning the atomic data in the initialization phase. CRC=EI electron impact collision, i.e., ionization, excitation or de-excitation or dissociation. (In older EIRENE version CRC=DS was also in use. This is now automatically identified with CRC=EI) CRC=CX (quasi-) resonant charge-exchange CRC=EL elastic collision CRC=PI in-elastic ion impact collision (not fully implemented) CRC=RC re-combination CRC=OT other MASSP Mass number (AMU) of the first particle involved in the collision, for the original data as being read from the data-file. EIRENE automatically re-scales data according to the mass number of the actual first particle involved in this particular collision process (e.g. if data are given for H atoms in the data-file, but are used for D or T atoms in an EIRENE run). MASST Mass number (AMU) of the second particle involved in the collision, for the original data as being read from the data-file. EIRENE automatically re-scales data according to the mass number of the actual second particle involved in this particular collision process. The next input deck is read only in case FILNAM = ADAS: ELNAME Name of element in ADF11 file, e.g. Fe for iron, C for carbon, W for tungsten, etc.... IZ Charge state in ADF11 file. This value is always between 1 and Zmax. For ionization rate coefficients IZ specifies the ionization rate coefficient for the ionization process from IZ−1 to IZ. For recombination process IZ specifies the recombination rate coefficient for the recombination process from IZ to IZ−1. The next input decks are read only in case FILNAM = PHOTON: IPRFTYPE IPLSC3 IMESS IFREMD NRJPRT to be written Then IFREMD cards are read II, KENN, IK6 II: labelling index, irrelevant KENN: to be written IR: to be written Note: The nomenclature used in reference [23] has lead to a certain confusion with regard to the masses of the particles involved in a particular collision process. A series of test calculations has revealed the following definitions implicit in the tables of references [23] and [24]: For the cross section energy scale, the relevant (laboratory frame) mass is that of the charged particle, assuming an ion- or electron beam incident onto a cold neutral gas at rest (i.e., the mass in the energy scale is neither the reduced mass nor is the energy given in units energy/AMU). For rate coefficients, depending upon both the beam energy and the target temperature, however, the neutral particles are considered to be the beam particles, i.e., their mass is the relevant mass for the energy scale now, whereas the charged particles are considered as target, with their mass being the relevant mass for the temperature scale. Following these conventions, EIRENE assumes that cross sections have been measured with the charged particle as projectile and the neutral target at rest. Thus MASSP is assumed to be the relevant mass number of the original data for the energy scale in cross sections (H.1) and for the temperature scale in the rate coefficients (H.2, H.3, H.4). If MASSP = 0 , MASSP is defaulted to the electron mass (AMU). I.e., MASSP is the mass number of those collision partners in the database, which would play the role of bulk particles in an EIRENE run, unless the automatic re-scaling to the true EIRENE bulk ions mass RMASSP(IPLS) of the bulk particle is carried out. MASST is assumed to be the relevant mass number for the projectile energy scale in H.3, H.6 and H.9, i.e., for the test particles in an EIRENE run. EIRENE then automatically converts the energy/temperature scales in the data fits to the correct scales for the particular particles masses (isotopes) involved in a collision event during the simulation process, if these do not happen to coincide with MASSP and MASST. DP (only for H.8, H.9, H.10) additional (e.g., potential) energy lost or gained, which is not already included in the rates. This may be needed due to the logarithmic fits used and if changes in sign arise. For example, in case of recombination of hydrogen ions, AMJUEL H.10, 2.1.8, DP=+13.6 must be specified. In reaction decks H.8, H.9 or H.10 in AMJUEL, this value of DP is now automatically read (and overwrites the value specified here). Hence: in the EIRENE versions later than June96 this input flag is redundant. The next flags control the options to extrapolate atomic data fits read from files. The default data in EIRENE and some reactions in the files AMJUEL already contain that information. EIRENE searches for the string ELABMIN and ELABMAX for cross section data. If it finds them, it uses them as RMN and RMX, respectively. In such cases the next card IFEXMN,.... and IFEXMX,.... may also be omitted, because the corresponding information is found from the atomic data file as well. The extrapolation expression is that of option IFEXM. = 5 (see below). Likewise, in case of reaction rates, EIRENE searches for TMIN, TMAX (for single parameter fits H.2, H.5, H.8 and H.11) and additionally for PMIN and PMAX (if H.3, H.4, H.6, H.7, H.9, H.10 or H.12) and the extrapolation parameters again are expected in the atomic data file. RMN minimum energy or temperature for data as read from file. Below RMN the data are extrapolated using the input card IFEXMN,FPARM(J),J=1,3. See below, and e.g., subroutine CROSS for the various options. Some cross section data in AMJUEL already contain the extrapolation information, and hence the default RNM=0 can be used there. RMX same as above, but for high parameter range extrapolation. IFEXMN, FPARM(J), J=1,3 various extrapolation expressions exlow(x) at the low end (x < RMN, x = Elab or x = T) are available: IFEXMN = 1 exlow(x) = 0 (e.g., cross section with nonzero threshold at RMN) IFEXMN = 2 y=x/FPARM(1) exlow(x) = FPARM(2) ·yFPARM(3) ·log(y) recommended for cross section extrapolation at high energy end, for reactions with nonzero threshold at FPARM(1) [eV]. IFEXMN = 3 y=log(x) ln(exlow(x)) = FPARM(1) + y ·FPARM(2) IFEXMN = 4 not in use IFEXMN = 5 y=log(x) ln(exlow(x)) = ∑I=13 FPARM(I) ·yI−1 IFEXMN < 0 linear extrapolation in log-log scale. Coefficients FPARM(J), J=1,2 are automatically determined, e.g., in function CROSS, FPARM(3) = 0.D0, and then extrapolation option IFEXMN = 5 is used with these coefficients (same as IFEXMN = 3 in this case). IFEXMX, FPARM(J), J=4,6 same as above, but for high parameter range extrapolation exhigh(x). *4A. Atoms Species Cards NATMI Total number of atomic species blocks IATM irrelevant; labelling index *4B. Molecules Species Cards NMOLI Total number of molecule species blocks IMOL irrelevant; labelling index *4C. Test Ions Species Cards NIONI Total number of test ion species blocks IION irrelevant; labelling index *4D. Photon Species Cards NPHOTI Total number of photon species blocks IPHOT irrelevant; labelling index The meaning of the variables in a species block is as follows: TEXT$
Name of the species on printout and plots
NMASS$Mass number NCHAR$
Nuclear charge number
NPRT$Flux units carried by one particle of this type and species. The "WEIGHT" of a test flight has dimensions: "atomic flux", i.e., "equivalent atoms per second". For example, a methane molecule CH4 should have NPRTCH4 = 5, if H atoms and C atoms have NPRTH = NPRTC = 1. NPRT is used to convert fluxes into equivalent atomic fluxes, for scaling, (see input flag FLUX in block 7), and for flux conservation in non-diagonal (species changing) events at surfaces and in the volume. (NPRT is irrelevant for atoms, and always set to one for them.) NCHRG$
Charge state (irrelevant for atoms and molecules, always set to zero for them)
[ISRF$] Species index flag for "fast particle reflection model" (see block 6A), irrelevant for molecules [ISRT$ ] Species index flag for "thermal particle re-emission model" (see block 6A)
ID1$not in use in versions 98 to 2002. Was "sputtered species index" in versions older than "98", but this species index has then become a surface property, see input block 3. In versions younger than 2002 this flag can be used to increase the number of secondary test particle groups from 2 or less (default) to 3, if ID1=3. This option then allows more complex fragmentation processes of large molecules as compared to the default. NRC$

> 0
total number of reaction decks to be read for this species
= 0
default model is to be used.
= -1
no collision processes for this particle species at all
NFOL$controls the motion of test particles. ≥ 0 default test particle tracing model is to be used. = -1 motion of test particle is not followed (static approximation). The test particle is destroyed immediately at it's point of birth by a collision. A collision estimator is used rather than a track-length estimator for all default volumetric tallies, as described in section 3.2. NGEN$
Maximum number of generations for this test particle species.
> 0
For each particle history the generation counter XGENER is reset to zero at birth, after surface events and after a non-diagonal event (modification of type and/or species of a particle in a collision event). It is increased at entropy producing events (elastic and diagonal charge exchange collisions). If XGENER > NGEN.. after a collision, no test particle secondaries are generated. The particle flux, parallel momentum flux and energy flux is put into the tallies PGENA, VGENA and EGENA for atoms, resp., PGENM, VGENM and EGENM for molecules, resp., PGENI, VGENI and EGENI for test ions, respectively and PGENPH, VGENPH and EGENPH for photons, respectively.
= 0
In case NGEN, (default) no generation limit is activated for atoms (NGENA=0), molecules (NGENM=0) or test ions (NGENI=0).
NHSTS$(only in versions 2004 and younger): NHSTS=-1 turns off the trajectory plot for this particular species. (default: NHSTS=0) ID3$
not in use
IREAC$Identification flag for collision data. IREAC$ = IR, and IR is the labelling index from the reaction card (see above).
The next three flags control number, species and type of particles involved in this particular collision process.
IBULK$pre collision bulk particle identifier = NLM M type flag ITYP for impacting bulk particle (electron (ITYP=5) or bulk ion (ITYP=4), irrelevant, defaulted to M = 4 for bulk ion impact collision and defaulted to M = 5 for electron impact collisions L irrelevant, defaulted to L = 1, because number of pre collision bulk particles is known from the type CRC of each collision process. N Species index of pre collision bulk ion (corresponds to IPLS in bulk ion species cards, see block 5). Irrelevant for electron impact collisions, CRC=EI. ISCD1$
first heavy secondary particle group identifier
= IJKLM
M
type flag ITYP for these first secondaries group (atoms (ITYP=1), molecules (ITYP=2), test ions (ITYP=3) or bulk ions (ITYP=4))
L
number of secondaries of this type and species
JK
Species index of secondary (IATM, IMOL, IION, or IPLS), two digits
I
relative importance of this first secondary group as compared to the second secondary group. If I ≤ 0, I is defaulted to I = 1
(currently not in use)
ISCD2$second heavy secondary particle group identifier = IJKLM. Same as ISCD1$, but for second secondary group.

Note:
The number of secondary electrons need not be specified, it is computed from charge conservation.
In case of collisions with only one heavy particle secondary (electron impact collision CRC=EI, or re-combination CRC=RC) ISCD1 is used and ISCD2 is irrelevant.
In case of charge exchange (CRC=CX) ISCD1 is used for the new state of the impacting bulk particle after the collision, whereas ISCD2 controls the options for the post collision state of that particle which was the test particle prior to the collision. Consistency of the particle masses of pre and post collision particles with this convention is checked in the initialization phase of each run.
In case of elastic collision (CRC=EL), the ISCD1 and ISCD2 flags are are irrelevant, because the colliding particles retain their identity by default.
In case of electron impact collision (CRC=EI), the ISCD1 and ISCD2 flags are symmetric, i.e. they can be interchanged without any effect on the calculation.
ISCDE$flags for post collision velocity space distributions for all particles involved in the collision event. Each of the five digits controls a different type of particle. Let us write JKLMN = ISCDE Let then furthermore the character # stand for J, K, L, M, or N respectively. Then J controls electrons, K controls pre collision (plasma) bulk ions (if any involved), L the secondary particles (if any). Each of these single integer flags controls the meaning of the energy parameters EP=EELEC (# [^=] J), EP=EBULK (# [^=] K), EP=ESCD (# [^=] L), for the electrons, for the pre collision bulk ions, and for the heavy secondary particles, respectively. # = 0 constant loss or gain of EP (eV) per collision event. # = 1 bulk temperature dependent loss -(1.5 ·T + EDrift) per collision, where T = Te for # [^=] J and T = Ti for # [^=] K; Te and Ti are the local electron and ion temperatures (eV) resp. at the point of collision, and EDrift is the kinetic energy of the drift motion (VXIN,VYIN,VZIN) of the background ions (only in case NLDRFT=TRUE), see input block 5. If # stands for a heavy particle secondary, i.e. # [^=] M and/or # [^=] N, then a fraction EP of the energy loss of the impacting bulk particle is distributed equally over all (if more than one) heavy particle secondaries of group 1 and 2 resp. # = 2 Let # stand for a pre collision bulk particle: # [^=] K or # [^=] L. The energy loss of the impacting bulk particle is calculated from the velocity and energy weighted rate coefficients. In this case EP must be the labelling index IR of the reaction card by which these moments have been read from the data-file. Let now # stand for a post collision heavy particle group, i.e # [^=] M and/or # [^=] N. Then the velocity vector of the secondary particle is sampled from a cross section weighted (and optional, if NLDRFT, shifted by the local drift velocity of the pre collision bulk particle) Maxwellian distribution at the local ion temperature Ti. The cross section is taken at the relative velocity of the impacting test particle and a particle sampled from the above mentioned (shifted) Maxwellian. This option permits to identify one particular collision partner from the background population, correctly accounting for the energy dependence of the cross section. # = 3 to be written # = 4 to be written # = 5 to be written IESTM Three digits IESTM=LMN. N, M and L are flags for the choice of estimators for particle, momentum and energy sink and source term contributions from the collision process. (Default: =0 track-length estimators.) IBGK Three digits IBGK=NML. Particle identifier for BGK self and cross collisions. Format as for IBULK, ISCD1 and ISCD2 species identifier flags. I.e., it has three digits, pointing to the type and the species index within that type-group. The second digit is irrelevant and defaulted to one (only binary collisions are considered). This flag controls options for non-linear test-particle - test-particle collision. Let ISPZ be the species index of a test-particle. If the collision is a self-collisions, this flag must point to the species ISPZ itself. Otherwise it must point to the 2nd test particle species ISPZ2, which is involved in the BGK-cross collision term. IBGK must be consistent with IBULK, which in this case is the "artificial background species" (= result from previous iteration) in the linearization of the collision operator. I.e.: if test particle species A collides with test particle species B (of same or of different species and type), then a bulk species ABbulk must be present in input block 5, such that species A and ABbulk have the same mass, nuclear charge and charge numbers. The density, flow field and temperature of collision partner ABbulk are iterated, using either the corresponding parameters estimated for species A (in case of self collisions) or using parameters obtained from the formulas for TAB, VAB (see section 1.9.1) for cross collisions between species A and B. The flag IBULK must point to the corresponding species IPLS=ABBULK in block 5. An example for input specifications of non-linear BGK collisions is given below in section 2.4.2. (Default: = 0, no iteration for non-linearity in collision kernel, no "artificial background species" involved in this collision, and no BGK-tallies are updated for iterating the artificial background species.) EELEC parameter EP for electron energy loss per collision EBULK parameter EP for pre collision bulk ion energy loss per collision ESCD1 parameter EP for velocity space distribution of secondaries EDUMMY (was earlier: ESCD2) in versions 1999 and older: parameter EP for velocity space distribution of second secondaries group. Now: this parameter is not in use anymore. Distribution of energy EP over secondaries is now by default inversely proportional to the secondary particle masses (from momentum conservation in collision). For backward compatibility of input files now: EP = ESCD1 + EDUMMY. FREAC multiplier to scale the cross section (and rate coefficients) of this process. Default: FREAC = 1.0. If FREAC ≤ 0, it is also defaulted to 1.0. To turn off this process, set FREAC to a very small number: e.g.: FREAC = 1.0000 E-30 FLDLM Flag for "generation limit" or "fluid-limit", Default: FLDLM = 0.0. This flag allows to cut trajectories after a certain maximum number of CX or EL collisions (generations) (FLDLM > 0.0) or if the mean free path compared to a geometric dimension falls below a certain value (FLDLM < 0.0) (fluid limit). The precise status of this option should be inquired from the code authors. 2.4.1 Default atomic and molecular data If default atomic data are requested for some test particle species, these are activated by setting the flag NRCM(IATM)=0 (or NRCM(IMOL)=0, NRCI(IION)=0, resp.) for a particular atom IATM (molecule IMOL or test ion IION), and by omitting the corresponding cards IREAC$,... and
EELEC$,... Then EIRENE tries to find, automatically the corresponding species types and indices of the collision products in the species blocks. It then tries to define all variables in the cards IREAC$,... and
EELEC,... for collision rates and collision kinetics, using default assumptions. Assume that the bulk ion H+ is specified as bulk ion species no. "a", D+ as bulk species "b", T+ as bulk species "c" and He+ as bulk ion species no. "d", then the default atomic collision models for H, H2, H2+ (and isotopes) and He are identical to a model that would result from the following specifications in input block 4. For atoms D (likewise for H or T), e.g., IATM=1, and atoms He (e.g., IATM=2): * ATOMIC REACTION CARDS NREACI= 5 1 HYDHEL H.2 2.1.5 EI 0 1 2 HYDHEL H.1 3.1.8 CX 1 1 3 HYDHEL H.2 2.3.9 EI 0 4 4 HYDHEL H.1 5.3.1 CX 4 4 5 HYDHEL H.1 6.3.1 CX 4 4 *NEUTRAL ATOMS SPECIES CARDS: NATMI SPECIES ARE CONSIDERED, NATMI= 2 1 D 2 1 1 0 1 -1 0 2 1 115 *14 0 00000 -1.3600E 01 0.0000E 00 2 *14 111 *14 00000 0.0000E 00 0.0000E 00 0.0000E 00 2 He 4 2 1 0 2 2 0 3 3 115 **14 0 00000 -2.4588E 01 0.0000E 00 4 **14 211 **14 00000 0.0000E 00 0.0000E 00 0.0000E 00 5 ***14 211 ***14 00000 0.0000E 00 0.0000E 00 0.0000E 00  Here ∗ is the species index of the D+ bulk ion in input block 5, ∗∗ is the species index of the He+ bulk ion in input block 5, and ∗∗∗ is the species index of the He++ bulk ion in input block 5. If a corresponding background species is not found in block 5, then the corresponding default process is de-activated. Hence: D and He atoms are ionized by electron impact, using the Corona data 2.1.5 and 2.3.9, respectively, from reference [23], and with constant electron energy loss per ionization of 13.6 and 24.588 eV, respectively. D atoms also undergo resonant charge exchange with D+ ions, utilizing the cross section data for reaction 3.1.8 (loc.cit.). The ISCDE parameters for CX processes are zero: = 00000 (or: = 00001, which has the same meaning because the fifth digit is presently irrelevant). The second digit (i.e., here = 0) controls the options for the velocity of the background ion going into the collision. For this choice the following model is activated: With respect to ions going into a CX collision the ion velocity distribution fi(vi) is assumed to be an isotropic, mono-energetic distribution , with vi0 = √{1.5 kTi}, then shifted by vDrift,i. The mean energy of ions undergoing charge exchange is taken to be 1.5·Ti eV plus the kinetic energy in their drift motion mi/2 (vDrift,i)2. This is used for the corresponding part in the CX momentum- and energy exchange rate coefficient. The reaction rate coefficient itself is evaluated as  〈σ·vrel〉 = σ(veff) ·veff with an effective relative velocity  veff = √ 3/2 kTi/mi + (v0 − vDrift,i)2 This approximation to the rate coefficient, for this particular choice of fi, is almost exact. It is assuming that the effect of the scalar product −vi·v0 in the relative velocity between ion- and neutral velocity becomes negligible after integrating σvrel over solid angle, which it would indeed exactly if σ ∝ 1/vrel. The velocity of the impacting ion is sampled, consequently, from a drifting, isotropic and mono-energetic distribution weighted by σ·|vrel|, with vrel the relative velocity between the two collision partners. Alternatively to the weighting mentioned above, a rejection method can be used to avoid the weighting and still simulate the same collision integral. Check in subr. VELOCX to find out which of the two physically equivalent methods are in use. Default reactions 5.3.1 and 6.3.1 (added in April 2006, to simplify input for pure He plasma simulations) are the two resonant charge exchange processes  He + He+ → He+ + He  He + He++ → He++ + He The weighting (or rejection-technique) was absent in versions younger than 99 for this particular type of (default) charge exchange collision integral approximation, which, strictly speaking, had therefore led to a slight violation of the second law of thermodynamics (H-Theorem), due to an inconsistency between the in-scattering term (determined by σCX and the out-scattering term (determined by 〈σvrel〉) in the Boltzmann CX collision integral, although mass and energy had been strictly conserved in each collision. Alternative to the choice ISCDE = 00000 made in the default model one may also use ISCDE = 01000 (with all the rest kept identical). In this case the rate coefficients are evaluated in the same way as described above, but neutrals emerging from CX collisions are sampled from the local drifting Maxwellian distribution of ion collision partners rather than from a drifting mono-energetic distribution. For hydrogenic molecules (H2, D2, HT,...) or molecular ions (H2+, D2+, HT+,....) default dissociation, ionization and dissociative recombination data are available. EIRENE uses the 6 reaction rates 2.2.5, 2.2.9, 2.2.10, 2.2.11, 2.2.12 and 2.2.14 from the data-set HYDHEL (loc.cit.) and again tries to identify the reaction products from the mass and nuclear charge number of the respective molecule IMOL (or test ion IION). In order to distinguish D2 from HT, is also uses the name TEXTM (or TEXTI), by looking for the appearance of the character "D" in the name (then: D2, or D2+) or for "H" or "T" (then: HT, or HT+). The default reaction kinetics are set as if they would have been specified by the following species blocks: * ATOMIC REACTION CARDS NREACI= 6 1 HYDHEL H.2 2.2.9 EI 0 2 2 HYDHEL H.2 2.2.5 EI 0 2 3 HYDHEL H.2 2.2.10 EI 0 2 4 HYDHEL H.2 2.2.12 EI 0 2 5 HYDHEL H.2 2.2.11 EI 0 2 6 HYDHEL H.2 2.2.14 EI 0 2 7 AMJUEL H.8 2.2.14 EI 0 2 . . . * NEUTRAL MOLECULES SPECIES CARDS: NMOLI SPECIES, NMOLI= 1 1 D2 4 2 2 0 0 1 0 3 1 115 113 0 00000 -1.5400E 01 0.0000E 00 0.0000E 00 0.0000E 00 2 115 121 000 00000 -1.0500E 01 0.0000E 00 3.0000E 00 3.0000E 00 3 115 111 *14 00000 -2.5000E 01 0.0000E 00 5.0000E 00 5.0000E 00 * TEST ION SPECIES CARDS: NIONI ION SPECIES, NIONI= 1 1 D2+ 4 2 2 1 0 1 0 3 -1 4 115 111 *14 00000 -1.0500E 01 0.0000E 00 4.3000E 00 4.3000E 00 5 115 *24 00000 -1.5500E 01 0.0000E 00 0.2500E 00 0.2500E 00 6 115 121 000 30300 7.0000E 00 0.0000E 00 7.0000E 00 0.0000E 00  As above, ∗ stands for the species index of the D+ bulk ion. Note that for reaction 6, the energy dependence of the cross section leads to a mean electron energy loss per event [ˉ(E)]el of about 0.88 ·Te, see reference [25]. This approximation is used in the default model for the electron energy loss Eel per collision and for the total kinetic energy release EK shared by the two product atoms. The reaction no. 7 above: AMJUEL H.8 2.1.14, provides values very close to this, by an independent fit to the correct electron energy loss for this process  〈σve Eel〉 = Eel ·〈σve〉 = kTe 〈σve〉· ⎛⎝ 3/2 + d ln〈σve〉 d ln(kTe) ⎞⎠ (4.2) The velocity distribution of the dissociation products is isotropic in the center of mass system (taken to be the system moving with the incident molecule), and the dissociation energy is shared by the dissociation products so as to preserve momentum. I.e., in case of unsymmetrical hydrogenic molecules the two dissociation products to not receive the same share of the dissociation kinetic energy release, but instead this is distributed inversely proportional to their masses. In case of mixed molecules, such as DT, HT or HD in case of default reaction models some rate coefficients are automatically splitt into two halfs and assigned to the proper product atomic particles. For hydrogenic bulk ions (H+, D+, T+) default volume recombination rates are available [26]. Further examples to be written 2.4.2 Neutral-Neutral collisions in BGK approximation In the example considered here we have the 4 neutral-neutral collision processes: D on D2, D on D and D2 on D2, D2 on D. In input block 4 there are three BGK relaxation rate coefficients (note: the rate coefficient for D2 on D is equal to the one for D on D2, hence 3 rate coefficients are sufficient).  . . some 20 other reaction cards in this example . 21 CONST H.2 EL 1 1 -2.1091E+01 0.2500E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 22 CONST H.2 EL 2 2 -2.0589E+01 0.2500E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 23 CONST H.2 EL 1 2 -2.0357E+01 0.2500E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00 0.0000E 00  In input block 4a (for atoms), we specify 4 collision processes for D atoms. Process no. 3 and 4 are neutral-neutral collisions. 1 D 2 1 1 0 1 -1 1 4 (i.e.: 4 processes to be specified for this D atom) . . 2 other process decks, e.g. CX and electron impact ionization . 21 2014 (=IBULK) 1001 0 111 (=IBGK , self collision) 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 23 2214 (=IBULK) 1001 0 112 (=IBGK, cross collision with IMOL=1) 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00  In input block 4b (for molecules), we specify 7 collision processes for D2. Process no. 6 and 7 are neutral-neutral collisions. 1 D2 4 2 2 0 0 1 0 7 (i.e.: 7 processes to be specified for this D2 molecule) . . . 5 other process decks. . . 22 2114 (=IBULK) 1001 0 112 (=IBGK, self collision) 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 23 2314 (=IBULK) 1001 0 111 (=IBGK, cross collision with IATM=1) 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00  In input block 5 there are at least 23 bulk ion species in this example. Of those species no. 20 to 23 are "artificial" BGK species for neutral-neutral self collisions specified above by IBULK = 2014, 2114, 2214 and 2314 in blocks 4a and 4b above.  . . . 19 other bulk ion species decks. . . 20 D(B) 2 1 1 0 1 1 0 0 21 D2(B) 4 2 2 0 1 1 0 0 22 DD2(B) 2 1 1 0 1 1 0 0 23 D2D(B) 4 2 2 0 1 1 0 0  2.4.3 Fitting expressions (IFTFLG) The default fitting expressions (IFTFLG=0) are single or double polynomial expressions, as described in [23] or [24] for the databases HYDHEL and METHANE, respectively. The following options are available in recent versions of the code: a) for cross sections, rate coefficients, rates: IFTFLG < 100 : rate coeff. or cross section, [cm3/s] or [cm2]. IFTFLG > 100 : rate = density * rate coeff, [s] (also used for spontaneous decay, e.g. radiative decay. mod(IFTFLG,100) = 10: only one fitting coefficient is read, i.e. the cross section, rate coefficient or rate is constant. b) for interaction potentials (elastic processes): to be written 2.5 Input for Plasma Background General remarks The background medium (mostly plasma) consists of NPLS (sometimes in the code: NPLSI) different so called "bulk particle" species, also referred to as "bulk ions", by abuse of language,. It is described by several blocks of NSBOX data each, (i.e. one datum per grid cell), namely • one array for the electron temperature Te (eV) (TEIN), • one common or NPLS (one for each bulk ion species) arrays for the ion temperature(s) Ti (eV) (TIIN), • NPLSI arrays (one for each bulk ion species), ion densities ni ( cm−3 ) (DIIN), • one common or NPLS (one for each bulk ion species) arrays for the (cartesian) drift velocities V = ( Vx , Vy , Vz ) (VXIN, VYIN, VZIN), either in cm/s or Mach number units. From the assumption of quasi-neutrality, the specified charge numbers of each ion species, and the data for ni, an array of NSBOX data for the electron density ne (cm−3 ) (DEIN) is computed (see next subsection: "derived background data"). The unit vector b = ( b1 , b2 , b3 ) (BXIN, BYIN, BZIN), parallel to the magnetic field is set to (0.,0.,1.) by default, i.e. pointing in z- (or toroidal-) direction. The default magnetic field strength is |B|=1 T in all grid cells. This default setting may be overwritten by using the INDPRO(5) flag and the corresponding input parameters (see below). Here either a (radial or x-) B-field "pitch" profile can be selected (INDPRO(5)=1,...4, and geometry level LEVGEO=1,2,3), with pitch = B2/Btot = b2. From this a cell centered B-field unit vector is internally constructed, by assuming that the first set of coordinate surfaces (x- or radial surfaces) are flux-surfaces (B-field aligned and constant pitch). The magnetic field direction then has components (0., pitch,√{1pitch2}) in the three coordinate directions. Currently only INDPRO(5)=3 profile option allows to define also the B-field strength (BFIN) [T], solely by the input flags described below. Alternatively and in all other options an external data source for the B-field vector in cartesian coordinates (INDPRO(5)=5,...,7, all geometry levels LEVGEO) is used directly. In this latter case the transfer of magnetic field data (BXIN,BYIN,BZIN and BFIN) for each cell into EIRENE proceeds from user supplied profile subroutine PROUSR (INDPRO(5)=5) or, from external file (INDPRO(5)=6,7), via work array RWK, and then from subroutine PROFR). The same applies to the "additional plasma profile array" ADIN, INDPRO(6), to permit usage of EIRENE output routines for plasma profiles such as pressure, energy fluxes, etc., which are not directly needed in an EIRENE run but which may be available from an external plasma code, here: INDPRO(6)=5,...,7 only. It also applies to the space and species dependent weight function WGHT (to be used for non-analog sampling techniques, INDPRO(7)), which is defaulted to 1.0D0 for all cells and all NSPZ species, and it applies to the zone volume array VOL cm−3, INDPRO(12). Defaults are the cell volumes computed in EIRENE subroutine VOLUME from the standard mesh data. Plasma data in the additional zones IADD, IADD=NSURF+1,NSURF+NADD outside the standard mesh are defaulted to the "EIRENE vacuum data" given below. Plasma data other than these in these zones have to be specified either explicitly in input block 8 (see below), or in the user supplied subroutine PLAUSR or with the INDPRO=7 option (see below, this section). Cell volumes in the additional cell region are defaulted to 1.0 cm−3, unless INDPRO(12)=4, 5, 6, 7. A standard mesh zone ICELL is automatically identified as a vacuum zone with regard to plasma species IPL (no collisions with the bulk particles IPL there, LGVAC(ICELL,IPL)=.TRUE., for electrons: IPL=NPLS+1) if at least one of the following conditions is met: for electrons: IPL=NPLS+1: TEIN(ICELL) ≤ TVAC , DEIN(ICELL) ≤ DVAC for bulk ion species IPL=IPLS: TIIN(ICELL,IPL) ≤ TVAC, DIIN(ICELL,IPL) ≤ DVAC and for all background species, IPL=0: LGVAC(ICELL,I)=TRUE for I=1,NPLSI and for I=NPLS+1. The EIRENE vacuum data are defaulted (in subroutine PLASMA) to: DVAC = 1.0000E 02 cm−3 TVAC = 2.0000E-02 eV Furthermore, zero plasma drift velocities are set in cells which are treated as vacuum zones. VVAC = 0.0000E 00 cm/s Note: A cell can be considered a vacuum cell with respect to electrons, (LGVAC(...,NPLS+1) = TRUE), without being a vacuum cell for all background ions. This is because, by abuse of language, also neutral particle species may be used as background "ion" species (NCHARP = 0), to include neutral-neutral collisions, e.g., by the iterative option NITER (input block 1). The Input Block *** 5. Data for Plasma background  NPLSI DO 51 J= 1, NPLSI  * read NPLSI species blocks with = P
 51 CONTINUE
(INDPRO(J), J=1,12)
IF (INDPRO(1).LE.5) THEN
TE0  TE1  TE2  TE3  TE4  TE5
ENDIF
IF (INDPRO(2).LE.5) THEN
(TI0(I)  TI1(I)  TI2(I)  TI3(I)  TI4(I)  TI5(I)  I=1,NPLSI)
ENDIF
IF (INDPRO(3).LE.5) THEN
(DI0(I)  DI1(I)  DI2(I)  DI3(I)  DI4(I)  DI5(I)  I=1,NPLSI)
ENDIF
IF (INDPRO(4).LE.5) THEN
(VX0(I)  VX1(I)  VX2(I)  VX3(I)  VX4(I)  VX5(I)  I=1,NPLSI)
(VY0(I)  VY1(I)  VY2(I)  VY3(I)  VY4(I)  VY5(I)  I=1,NPLSI)
(VZ0(I)  VZ1(I)  VZ2(I)  VZ3(I)  VZ4(I)  VZ5(I)  I=1,NPLSI)
ENDIF
IF (INDPRO(5).LE.5) THEN
B0  B1  B2  B3  B4  B5
ENDIF
IF (INDPRO(12).LE.5) THEN
VL0  VL1  VL2  VL3  VL4  VL5
currently option INDPRO=6 and INDPRO=7 are not available for volumes.
ENDIF



Meaning of the Input Variables for Plasma Parameters

The meaning of the variables in the bulk ion species cards is as in block 2.4 for the species cards specified there.
New options since 2003:
However, after the ID3 flag, there may be 2 additional input flags: (only for Eirene2003 and younger):
..... CDENMODEL NRE


Format: ....,1X,A10,1X,I2.
These flags, if included, control additional options to set the density, temperature and flow field of the selected species IPLS. If the string CDENMODEL is found here for a species IPLS, then, after reading the species and reaction cards for IPLS, one (default) or NRE further input cards are expected.
CDENMODEL = 'Multiply'
docu to be written
CDENMODEL = 'Constant'
docu to be written
CDENMODEL = 'fort.13'
Read data for this species from file fort.13.
ISP, ITP


ISP is the number of a bulk species on the file fort.13 (e.g. written in an earlier EIRENE run, see NFILEL option in input block 1). The parameters (fields of density, temperature, flow velocity) for species IPLS are set from those of species ISP on fort.13
ITP is the type of the particle, i.e., always: ITP=4 for this option. Internally set: ITP=4, bulk particle type, independent of input value for ITP.
CDENMODEL = 'fort.10'
Read data for this species from file fort.10.
ISP, ITP, ISTR


ISP, ITP is the number and type of a test particle species on the file fort.10, respectively, from stratum no. ISTR (i.e. written onto fort.10 either in the present run or an earlier EIRENE run, see NFILEN option in input block 1). The parameters (fields of density, temperature, flow velocity) for species IPLS are derived from from those of species ISP, ITP, ISTR on fort.10.
Note: at the end of a full EIRENE run the corresponding bulk particle profiles are reset to those evaluated in this run, for further use in post processing routines such as the diagnostic module (see Section 2.11), printout, plotting or iterative mode (input flags NITERI, NITERE, and user subroutine MODUSR, see Section 2.1)
CDENMODEL = 'Saha'
Set parameters from Saha-equilibrium of ionization states. The additional card (only one) reads:
ISP, ITP, ISTR, ....


to be written
CDENMODEL = 'Boltzmann'
Set parameters from Boltzmann-equilibrium of excited states. The additional card (only one) reads (format: 3I6,6x,2E12.4)
ISP, ITP, ISTR, G_BOLTZMANN, DELTA_E


docu to be written
CDENMODEL = 'Planck'
(only for background photons) to be written
CDENMODEL = 'Corona'
Set parameters from Corona-equilibrium. The additional card (only one) reads (format: 3I6,1X,A6,1X,A4,A9,A3,E12.4):
ISP, ITP, ISTR, FILNAM , H123 , REACTION, CR, A_CORONA


ISP, ITP is the number and type of a particle species, either one of the already defined bulk species IPLS: ISP=IPLS, (for ITP=4), or from the file fort.10 (for ITP=0,1,2 or 3) from stratum no. ISTR (i.e. written onto fort.10 either in the present run or an earlier EIRENE run, see NFILEN option in input block 1). The parameters (fields of density, temperature, flow velocity) for species IPLS are set from those of species ISP, ITP, ISTR.
This model then constructs a background species IPLS which is in corona-equilibrium with these ISP,ITP, ISTR particles.
set parameters from collisional radiative equilibrium, using reduced population coefficients from atomic database. The NRE additional cards read:
ISP, ITP, ISTR, FILNAM  H123 ....


ISP, ITP is the number and type of a particle species, either one of the already defined bulk species IPLS: ISP=IPLS, (for ITP=4), or from the file fort.10 (for ITP=0,1,2 or 3) from stratum no. ISTR (i.e. written onto fort.10 either in the present run or an earlier EIRENE run, see NFILEN option in input block 1). The parameters (fields of density, temperature, flow velocity) for species IPLS are set from those of species ISP, ITP, ISTR.
This model then constructs a background species IPLS which is the sum of NRE species in collision-radiative-equilibrium with these ISP,ITP, ISTR particles.
INDPRO
Flag-array for the type of profile. The last digit (between 1 and 9) controls the type of profile. A second digit and/or the sign control further options, as described below. INDPRO is an array of length 12. Each element in this array controls one particular input tally, namely:
INDPRO(1)
for TEIN
INDPRO(2)
for TIIN By the default (0 < INDPRO(2) < 10): one common Ti profile for all "bulk ion" species is set. I.e., read only one profile card.
New options since 2001:
Values of INDPRO(2) larger than 10: use only the last digit, and one Ti profile card must be read for each bulk ion species IPLS. For example: INDPRO(2)=15 or =25, means: one separate ion temperature must be specified for each bulk ion species, and the profile type is 5. (INDPRO(2)=-5 would do the same.)
INDPRO(3)
read NPLSI cards, for DIIN , IPLS = 1,NPLSI
INDPRO(4)
read NPLSI cards, for VXIN, VYIN, VZIN , IPLS = 1,NPLSI
New options since 2001:
By the default (0 < INDPRO(4) < 10): one separate flow field for each bulk species is set, velocity is given in cm/sec.
Negative value of INDPRO(4) means: Mach number units instead. The sound speed cs is taken to be the isothermal ion acoustic speed of species IPLS.
ABS(INDPRO(4) is then used as flag for the choice of profile type.
Values of INDPRO(4) larger than 10: use only the last digit, and only one common flow field is set for all bulk ion species. Note the difference to the Ti options: there the meaning of INDPRO larger than 10 was exactly opposite to the meaning here for the flow fields (due to historical reasons and for backward compatibility of input files. EIRENE had originally by default one single common ion temperature, but one flow field for each bulk ion species).
INDPRO(5)
for PITCH (later, in initialization phase, converted into cartesian unit B-field vector BXIN,BYIN,BZIN)
Pitch is defined as By/Btot (LEVGEO=1), Bθ/Btot (LEVGEO=2), or as Bpol/Btot (LEVGEO=3), where Bpol is the direction along the polygons.
New options since 2001: In case INDPRO(5)=3 (flat profile) the two redundant input parameters B2, B3 are used to define a constant B-field strength [T], see below under profile type INDPRO=3.
INDPRO(6)
INDPRO(7)
for WGHT (to be written)
INDPRO(12)
for VOL.
For the profiles no. 6 and 7 only the options INDPRO=5 or INDPRO=6 exist.
For the profile no. 12 (cell volumes) the options INDPRO(12) = 4, 5, 6, or INDPRO(12) = 7 are active options. For all other values of INDPRO for these latter four profiles the default profiles described above ("General remarks") are set.
For each profile up to 6 parameters P0 , ..., P5 are read, e.g. TE0, ... , TE5 for the Te -profile, TI0, ... , TI5 for the Ti -profile, and so on.
Depending upon the value of INDPRO one of the profile routines PROFN, PROFE, PROFS, ... etc. is called from subroutine PLASMA.
INDPRO = 1-4

NR1STM plasma data are defined on the one- dimensional zone-centered grid
RHOZNE(J), J=1, ...,NR1STM (NR1STM = NR1ST - 1 ) .
Vacuum data are specified in zone NR1ST.
These NR1ST data are copied NBLCKS times to define
NBLCKS identical profiles, totally : NBLCKS · NR1ST = NSURF data (subroutine MULTI).
Here the number of copies is NBLCKS = NMULT · NP2ND · NT3RD. see "Input Data for Standard Mesh" (block 2).
INDPRO = 1
(see subroutine PROFN)
RHOSRF(1) ≤ x ≤ P5:
P(x) = P1+( P0 − P1 ) · ( 1 −( [(xRHOSRF(1))/(P5 − RHOSRF(1))] ) P2) P3

in particular :
 x = RHOSRF(1) → P(x) = P0
 x = P5 → P(x) = P1

P5 ≤ x ≤ RHOSRF(NR1ST) :
P(x) = P1 · exp(- (x - P5) / P4)
INDPRO = 2

(see subroutine PROFE)
RHOSRF(1) ≤ x ≤ P5 :
P(x) = P0 ·exp(- (x - P1) / P2)
in particular : x = P1 → P(x) = P0
P5 ≤ x ≤ RHOSRF(NR1ST) :
P(x) = P(P5) ·exp(- (x - P5) / P4)
(P3 : no meaning)
INDPRO = 3

(see subroutine PROFS, PROFS(P0,P1,P5,PVAC)
RHOSRF(1) ≤ x ≤ P5 : P(x) = P0
P5 ≤ x ≤ RHOSRF(NR1ST) : P(x) = P1
Parameters P2, P3, P4 have no meaning, except for the magnetic-field (INDPRO(5)=3): Here PROFS(P0,P1,P5,PVAC) is used to define the pitch angle profile, and it is also used to set the B-field strength (Tesla) by calling PROFS(P2,P3,P5,PVAC)
INDPRO = 4

INDPRO = 5

EIRENE calls a user supplied plasma profile routine:
Subroutine PROUSR (RHO, INDEX, P0, P1, P2, P3, P4, P5, PVAC, NDAT).
NDAT data must be defined on
RHO(J), J=1,NDAT e.g. by using the profile parameters P0,...,P5 and any user supplied profile function. See section 3.6.
(NDAT=NSURF)
The flag INDEX is as for the "INDPRO = 6" option, see below, and sections 3.6 and 4 . It determines, which of the background medium profiles has to be provided at a particular call to subroutine PROUSR.
PVAC is the EIRENE vacuum value for the profile in question.
INDPRO = 6

EIRENE calls a the routine:
Subroutine PROFR (RHO,INDEX,N1I,N1,NDAT).
N1I (e.g., NPLSI) profiles, of NDAT= NSURF data each, are read onto RHO(I1I,J), J=1,NDAT, I1I=1,N1I from a work array RWK. N1 (e.g., NPLS) is the leading dimension of the array PRO as specified in the calling program. The data must be written onto array RWK in the initialization phase, e.g., in subroutine INFCOP or MODUSR. The location of a particular profile on this work array is determined by the value of INDEX, see section 4 below. By this option plasma parameters may directly be transferred into EIRENE from other files, e.g. from plasma transport codes.
NDAT (= NSURF) is the number of cells in the "Standard Mesh". Parameters in the additional cell region
are set to the "EIRENE vacuum values"
For further details of subroutine INFCOP see 2.1 "Input Data for operating mode", input variable NMODE and sections 2.14 and 4 for an example.
INDPRO = 7

Same as INDPRO=6 option, except that NDAT=NSBOX rather than NDAT=NSURF. Hence in this case background data are read from RWK also for the additional cell region

2.5.1  Derived Background Data

Given the data describing the background medium (plasma) in each cell (on common block COMUSR) a set of further, "derived background data" profiles is computed and stored also on COMUSR. This is done by a call to subroutine PLASMA_DERIV. The following quantities are defined there:
DEIN
Electron density (cm3), derived from all NPLSI ion densities and the charge states NCHRGP(IPLS) of each ion species.
EDRIFT
Kinetic energy (eV) in drift motion, for each background species IPLS, derived from the mass MASSP(IPLS) and the flow velocity VXIN,VYIN,VZIN.
LGVAC
Vacuum flag (logical) to identify regions, in which no collision data for all or certain background species are computed, i.e., the collisionalities are set to zero in these cells for those background species. See also explanation above, with respect to the EIRENE vacuum data TVAC,DVAC,VVAC.
NSTGRD
flag for special treatment of some selected cells
=0
default
=1
This cell is a dead cell, not to be seen by the particles (isolated from the computational domain. NSTGRD(ICELL)=1 can be specified in problem specific routines (...USR) or in routines interfacing EIRENE with external codes: INFCOP (code segment: COUPLE....)
=2
indicates that this cell belongs to a grid cut, i.e., is not a valid cell for particle tracing.
=3
indicates that this cell is not a real cell but only the storage position for spatially averaged tallies. E.g.: cell no. NR1ST, 2*NR1ST,....etc. .

2.6  Input Data for Surface Interaction Models

General remarks

An outline of surface interaction models in general terms was given in section 1.4. As for the implementation of such models in EIRENE, there are two parts to the surface interaction data block. The first part contains data which are general to the EIRENE reflection model ("Block for General Reflection Data"). The second part may be different for each surface element and thus must be specified for each reflecting "non-default" surface of the "standard mesh" and for each reflecting "additional surface". It is referred to as "Block for Local Reflection Data". If this block is missing and if the surface MSURF is neither transparent (ILIIN < 0) nor purely absorbing (ILIIN=2) nor a "mirror surface" (ILIIN=3) nor a "periodicity surface", (ILIIN ≥ 4) then the default reflection model (FeTarget, ... ) is activated for this surface.
The surface interaction model comprises the following information:
1. pf (RPROBF): probability for reflection of a "fast" particle
(= 0. for incident molecules),
2. pt (RPROBT): probability for re-emission as "thermal particle"
(with surface temperature, see the flags EWALL)
3. pa (RPROBA): probability for absorption
(surface sticking, pumping, etc...)
4. ps : probability for (physical) sputtering
5. pc : probability for (chemical) sputtering
6. For the reflected "fast particles" and for the sputtered particles the type, species and the distribution of energy, polar and azimuthal angle must be given.
7. For the re-emitted thermal particles the type, species and only the distribution of energy must be given; the angular distribution follows the cosine law by default.
All these data may be functions of incident energy Ein, of incident polar angle θin against the surface normal and of projectile and target species. The three probabilities pf, pt and pa must add up to 1; ps and pc are sputtering yields per incident particle and, thus, may occasionally be larger than 1. They are called a probability here only by abuse of language.
All this information is contained in each of three so called "Reflection Models", namely the "Database Reflection model" (see references [15], [16]), the "Modified Behrisch Matrix model" (references [13], [27]) or the "User Supplied Reflection model" (see section 3.3), and in several so called "Sputter Models" ([28], [29], [30].
These models may, to some extend, be modified by the input flags described here. For example, a surface may be purely absorbing (for all test particles incident onto it, ILIIN = 2 option), acting like a mirror for the neutral test particles, by the ILIIN = 3 option (see section 2.3A and 2.3B), or enforce periodicity, e.g., because of symmetry conditions (ILIIN ≥ 4 option). In these two latter cases all test particles of all species are reflected with the species unchanged, with probability one and elastically: Eout = Ein. In case of the "mirror reflection option" the normal component vn of the incident velocity vector v is reversed: vn → −vn. In case of a "periodicity surface" both the position and the speed unit vector are altered according to the specific periodicity of the configuration.
If, instead, a surface reflection model is chosen (ILIIN = 1), then some of the input flags described here (RECYCF, RECYCT, RECYCS, RECYCC) as well as the flags (RINTEG, EINTEG, AINTEG) may be used to modify conveniently the otherwise rather unhandy reflection data and formulas.

The Input Block for General Reflection Data

*** 6A. Data for General Reflection Models
 NLTRIM
Path-Card" for TRIM surface reflection files"
(machine specific, may be omitted)


If a "Path Card" is included, then read an arbitrary number of "Target-Projectile Specification Cards" of the format
 A_on_B


Up to NHD6 (see PARAMETER statements, section 3.1) data files for different Target-Projectile combinations can be read.
If no "Path Card" has been specified, or after the "Target-Projectile Specification Cards", continue reading:
 (DATD(I) I=1,NATM)
(DMLD(I) I=1,NMOL)
(DIOD(I) I=1,NION)
(DPLD(I) I=1,NPLS)
ERMIN  ERCUT  RPROB0  RINTEG EINTEG AINTEG


optionally, an arbitrary number of surface models, labelled by the character string ′modname′, may be read next.
SURFMOD_modname
ILREF  ILSPT ISRS  ISRC
ZNML EWALL EWBIN TRANSP(1,N) TRANSP(2,N) FSHEAT
RECYCF  RECYCT  RECPRM  EXPPL  EXPEL  EXPIL
RECYCS  RECYCC  SPTPRM  (this line may be omitted,
then: default sputter model, see below)
next: within each SURFMOD-sub-block there may be an
arbitrary number of lines
VARNAME  SPZNAME   VALUE
(format: CHARACTER*8,X,CHARACTER*8,*)



Meaning of the Input Variables

NLTRIM
TRIM database is used, if "Database Reflection Model" is specified in at least one block for local reflection data. Data are read from data-set FT21 (no "Path Card" specified, old option), or from the domain specified by the "Path Card" (new option, see next card). If the old option is used, then the complete TRIM file is read, containing the first 12 TRIM target-projectile combination data-sets listed in section 1.4, i.e., the files H_on_Fe to T_on_W.
If the parameter NHD6 < 12, (section 3.1) then only the first NHD6 files are read from FT21.
"Path Card"
This card is machine specific. If EIRENE finds a card containing the string 'PATH' or 'path', it assumes that this card specifies the path to the domain containing the TRIM surface reflection data files.
A_on_B
Name of a particular TRIM data file in the domain specified by the path card. E.g., H_on_Fe would include the data file for hydrogen onto iron into the EIRENE run. Up to NHD6 such "Target-Projectile Specification Cards" may be included. For a complete list of such files currently available see again section 1.4.
Note: if during a Monte Carlo simulation a projectile A hits a target B, for which no TRIM data file has been specified, but still NLTRIM=TRUE, then EIRENE searches the "closest" of all its TRIM data files (with respect to a reduced mass argument) and uses this target-projectile combination together with a reduced mass scaling of incident particle energy. Hence, a TRIM data set is chosen such that the reduced mass scaling factor is as close to one as possible amongst the files available.
DATD
distribution for sampling the species index of reflected or otherwise emitted atoms. The NATMI relative frequencies
DATD(IATM) IATM = 1,NATMI
are used to produce the corresponding cumulative distribution DATM in order to facilitate sampling (inversion method).
DMLD
as for DATD, but for molecules. Cumulative distribution is DMOL.
DIOD
as for DATD, but for test ions. Cumulative distribution is DION.
DPLD
as for DATD, but for bulk ions. Cumulative distribution is DPLS.
The activation of these distributions at surface events during the particle history generation process is controlled by the surface species index flags ISRF$, ISRT$ read in the species sub-blocks of section 2.4 and 2.5. By default the following conventions are used (and can be overruled by calls to subroutines REFUSR, SPTUSR only, see section 3.3):
ISRF$fast particle reflection species flag: > 0 If ISRF$ ≤ NATMI, then ISRF$=IATM, the species labelling index for the reflected fast atom. If ISRF$ > NATMI not in use, warning and error exit.
= 0
no fast particle reflection for this species: pf = 0
< 0
not in use, warning and error exit.
ISRT$thermal particle re-emission species flag: incident atoms, test ions and bulk ions: > 0 If ISRT$ ≤ NATMI, then ISRT$= IATM, the species labelling index for the re-emitted thermal atom. If ISRT$ > NATMI, then IATM is sampled from the distribution DATM(IATM).
= 0
no thermal particle re-emission for this species: pt = 0
< 0
molecule is re-emitted, and -ISRT$is used to identify the species labelling index IMOL for the re-emitted molecule, exactly as described above for re-emitted atoms. The relevant sampling distribution for species index IMOL in case -ISRT$ > NMOLI is DMOL(IMOL).
incident molecules:
> 0
If ISRT$≤ NMOLI, then ISRT$ = IMOL, the species labelling index for the re-emitted thermal molecule.
If ISRT$> NMOLI, then IMOL is sampled from the distribution DMOL(IMOL). = 0 no thermal particle re-emission for this species: pt = 0 < 0 not in use, warning and error exit These latter two surface-species index flags are considered to be particle properties, hence they are read in the particle specification blocks 2.4 and 2.5. The flags in the next card can be used to modify the pre-programmed fast particle reflection models, to some extend at least. ERMIN For incident particle energies below ERMIN, the "fast" particle reflection model is switched off. Only the "thermal" particle model is used. ERCUT,RPROBF These variables may be used to modify the default "Behrisch Matrix" reflection coefficients for particles incident on a surface at low energies Ein. The original data [13] are used only for Ein > ERCUT and for normal incidence ϑin = 0 . In the range ERMIN < Ein < ERCUT the particle reflection coefficient pf (Einin = 0) is replaced by a smooth cubic interpolation curve pf (Ein) such that pf (0) = RPROBF . The original "Behrisch Matrix" is recovered by setting ERCUT ≤ 0. The next three flags modify the particle, energy and momentum reflection coefficients, and/or the resulting distributions in post-reflection energy and angle, respectively. These flags only apply to incident particles of type 1, 3, 4 and 5, i.e., not for incident molecules because for those RPROBF = 0 by default. RINTEG > 0 Fixed (independent of energy and angle of incidence) particle reflection coefficient. The fast particle reflection probabilities RPROBF are set to pf=MIN(1−pa,RINTEG), regardless of the reflection model selected by the flag ILREF in block 6B. pa is kept as specified, and pt is then recomputed as pt = 1 − pfpa. RINTEG ≥ 1 − pa enforces the fast particle reflection model for all un-pumped incident particles, i.e., RINTEG is internally reset to 1 − pa. = 0 Default: fast particle reflection probability as defined by the reflection model chosen. < 0 The fast particle reflection probabilities RPROBF are set to 1.0, i.e., even pumping is turned off (as distinct from the choice RINTEG=1.0, which would preserve the pumping speed at a surface). EINTEG > 0 Fixed (independent of energy and angle of incidence) energy reflection coefficient. This choice replaces the energy random sampling procedure in the fast particle reflection model by an reflection assumption: Eout = Ein ·EINTEG. = 0 Default: no modification of energy distribution in the reflection model for fast particle reflection. < 0 elastic (Eout = Ein) reflection for all particles reflected according to the fast particle reflection model. Hence: same as EINTEG=1.0 AINTEG (not ready to use) > 0 Fixed (independent of energy and angle of incidence) momentum reflection coefficient. This choice replaces the angle random sampling procedure in the fast particle reflection model by a momentum reflection assumption such that, on average,: vout = vin ·AINTEG. = 0 Default: no modification of angular distributions in reflection model for fast particle reflection. < 0 not in use Note: Some "General Reflection Model" data of input block 6A are copied identically NLIMI + NSTSI times onto appropriately dimensioned arrays (with the same names) in order to make them dependent upon the surface labelling index as well. The first index is then the species index, the second index is the surface number. Presently this "localization option" is available for the flags to be written  RF0USR or SP0USR of the user supplied reflection routine REFUSR; RF0USR and SP0USR are called in the initialization phase of an EIRENE run (see section 3.3). By a similar strategy, some of the species independent flags in the next sub-block 6B can be made species-dependent. See, for example, RECYCF, RECYCT, RECYCS and RECYCC below. *** 6B. Data for Local Reflection and Sputtering Models The next 3 (or 4) lines comprise the sub-block for local reflection data. Such sub-blocks can be included in each "surface deck" (block 3A, 3B) to overwrite the default reflection model for each individual surface. ILREF ILSPT ISRS ISRC ZNML EWALL EWBIN TRANSP(1) TRANSP(2) FSHEAT RECYCF RECYCT RECPRM EXPPL EXPEL EXPIL RECYCS RECYCC SPTPRM (this line may be omitted then: default sputter model, see below) An arbitrary number of such sub-blocks of 3 (or 4) lines each may also be defined under a certain label "SURFMOD_modname" and be included in the input file in block 6 directly after the ERMIN,... deck In blocks 3a and 3b surfaces may be assigned a particular local reflection model by a card reading SURFMOD_modname. Many surfaces may then be linked to the same surface reflection sub-block. Example: SURFMOD_BERYL_SPT_300K 1 2 0 0 9.04000E+02-2.60000E-02 0.00000E+00 0.00000E+00 0.00000E+00 2.80000E+00 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00 5.00000E-01 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00 SURFMOD_CARB_SPT_1153K 1 12 0 0 1.20600E+03-1.00000E-01 0.00000E+00 0.00000E+00 0.00000E+00 2.80000E+00 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00 5.00000E-01 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00 TRANSP1 D2 0.5000E 00 TRANSP2 T2 0.2500E 00 RECYCF D 0.0000E 00 RECYCT D2 0.9500E 00 SURFMOD_CARB_SPT_812K 1 2 0 0 1.20600E+03-7.00000E-02 0.00000E+00 0.00000E+00 0.00000E+00 2.80000E+00 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00 5.00000E-01 1.00000E+00 1.00000E+00 1.00000E+00 1.00000E+00  Meaning of the Input Variables ILREF Flag for choice of local reflection model = 1 TRIM database reflection model is used. NLTRIM must be .TRUE.. = 2 "modified Behrisch Matrix model" is used = 3 user supplied reflection model (see section 3.3: Subroutine REFUSR) Default: ILREF = 2 ILSPT Flag for choice of local sputtering model Let ILSPT = MN, with M and N single digit integers each. Then N controls the options for physical sputtering, and M controls chemical sputtering. See subroutine SPUTER. N = 0 no physical sputtering at this surface N = 1 constant physical sputtering rate (see parameter RECYCS below) N = 2 modified Roth-Bogdansky formula for sputter yield, Thompson energy distribution and cosine angular distribution for emitted particles (see references [29] and [30]). N = 9 (was option N=3 in Eirene-2004 and older) user supplied sputtering model (see section 3.3: Entry SPTUSR to subroutine REFUSR) M = 0 no chemical sputtering at this surface M = 1 constant chemical sputtering rate (see parameter RECYCC below) M = 2 "Roth formula" for chemical sputter yield, thermal distribution for emitted particles (see reference [31]), "weak flux dependence option A6". M = 3 "Roth formula" for chemical sputter yield, thermal distribution for emitted particles (see reference [31]), "strong flux dependence option A7". M = 4 "Roth formula" for chemical sputter yield, thermal distribution for emitted particles (see reference [31]), "new flux dependence option A8 (2004)". M = 5 not in use M = 6 "Haasz-Davis 1998 formula" for chemical sputter yield M = 7 "Haasz-Davis 1998 formula" for chemical sputter yield, and multiplicative factor for flux dependence (Roth, 2004). M = 9 (was option N=3 in Eirene-2004 and older) user supplied sputtering model (see section 3.3: Entry SPTUSR to subroutine REFUSR) Default: ILSPT=0 The next two surface-species index flags ISRS$ and ISRC$control the species of sputtered particles. Hence they are considered surface properties and are read in the surface specification decks in blocks 2.3.1, 2.3.2 and/or 2.6. ISRS$
sputtered particle species flag (physical sputtering).
> 0
both the sputtered particle and the reflected particle (if any) will be followed. Their contribution to surface particle and energy fluxes is stored in surface averaged tallies 1 to 24, i.e., sputtered particles are not explicitly distinguished from reflected particles in the balances. Furthermore the "sputtered flux surface tallies" 25 to 28 are updated. The species index of the sputtered particle (atom) is IATM=ISRS. Hence, on input, ISRS ≤ NATMI.
= 0
There is no physical sputtering for the particular surface element and incident species (note: ISRS=ISRS(ISPZ,MSURF), i.e., ps = 0 here.
Only the reflected particles are followed. Their contribution to surface particle and energy fluxes is stored in surface averaged tallies 1 to 24.
< 0
Same as ISRS > 0, however, the species index of the sputtered particle is determined automatically from comparing the charge and mass numbers of the atomic species (input block 4a) with the corresponding data of the surface element.

3.2.3  Snapshot estimated volume averaged tallies, UPNUSR

By default snapshot tallies are only scored in time dependent mode, see input flags in Section 2.13. They can, however, also be used as a third alternative to tracklength and collision estimators for any stationary response. In particular for diffusive processes, as e.g. ion transport with Fokker Planck type Coulomb collision operators, such "snapshot averaging" is the most frequent option used in Monte Carlo applications to estimate stationary moments. See paragraph B below.

A: time dependent estimates

The default snapshot estimated tallies at a fixed time t*, see Section 2.13, are updated in subroutine TIMCOL. Currently these are only the census tallies and the particle- and energy fluxes at census. Subroutine TIMCOL is called from the particle tracing subroutines FOLNEUT and FOLION (defaults), in case of time dependent mode (see input block 13). Non default snapshot tallies can be scored in UPNUSR (additional user supplied tallies, called from TIMCOL). Here we have
 WSNAP = k ^ wi .
In EIRENE variables:
 WSNAP = WEIGHT     .
Then, in order to estimate a response for a detector function g on snapshot tally number ISNV, a statement in subroutine UPNUSR should read:
 SNAPV(ISNV,ICELL) = SNAPV(ISNV,ICELL) + WSNAP * g(...)  .


E.g., with g=1 this is an unbiased estimate of cell averaged particle density at time t=t*, see (3.2), i.e., just counting the particle weights at t=t* provides an unbiased estimate of the particle density. This is intuitively clear, but also mathematically correct since these snapshot estimators are formally derived as special cases of tracklength- collision or surface flux estimators (interpreting the time-horizon sub-manifold as "time-surface", by abuse of language) by including in g a delta function in time : g(t,r,v)=δ(t*)[(g)\tilde](r,v).
Subroutine UPNUSR must at least include module COMPRT (for the particle weight WEIGHT). Here g is evaluated at the point of time t* at which TIMCOL is called. Scaling of snapshot tallies is done by interpreting the linear source strength "FLUX" scaling factor (given in input block 2.7 for each stratum) as a total number of particles, rather than as flux (particles/s), i.e. "FLUX is the integral of the source distribution not only over physical space, but also over time interval DTIMV from initial time t0 to time t* (Section 2.13).

B: stationary snapshot tallies

A stationary value of snapshot tallies, to be compared with other tallies in stationary mode, is obtained by setting the flag NTMSTP of block 2.13 to negative values. Stationary results in this "time dependent mode" are obtained by launching all test flights at time t0, but scoring snapshot tallies not only at t*=t1=t0+DTIMV, but also at t2=t0+2 DTIMV, t3=t0+3 DTIMV, ..., i.e. at all times tn=t0+n DTIMV until the flight terminates. If the source distribution and the background medium parameters are constant in time, then the snapshot score contribution from time ti can be interpreted as score at t*=t1 resulting from the source at earlier time t0−(i−1) DTIMV. Hence this stationary snapshot score at t* is in fact an integration over all contributions from earlier times t < t0.
When scaling the time- snapshot tallies at t=t* in this stationary mode to absolute units, then an extra multiplicative factor DTIMV (i.e., the time interval between subsequent scores along a particle history) is applied for snapshot tallies, because the stationary (time-averaged) volume averaged tallies are scaled interpreting the stratum source strength FLUX as flux: "atomic particles per second", whereas the snapshot tallies require this same scaling factor FLUX to be the absolute number of (atomic) particles (flux integrated over one time interval from t0 to t1=t0+DTIMV).
Further considerations regarding the scaling factors for volume averaged tallies in EIRENE, with respect to dimensionality of phase space of the problem (1D, 2D, 3D, stationary or time-dependent) are given in Section 1.3.2.

3.2.4  Surface averaged tallies, UPSUSR

Finally we note that "surface averaged tallies" (table 5.3) may be considered just as special cases of the "volume averaged tallies" described above, if appropriate use is made of δ -functions. Let therefore ρ be a co-ordinate normal to a surface S at the strike point rS of a test flight from ri to ri+1 with speed unit vector Ω = v /(| v |). I.e., the surface is described locally by the equation ρ = 0 at the point rSS, at which we may define an orthonormal basis eρ, eρ′ ,eρ". Furthermore, rS = ri+1 , which means that one may treat surface events exactly as collisions with the background medium in our terminologies. Then a "surface response function" gS defined as

 gS (j, r, v) = δ(ρ) ·g(j, r, v)
is the detector function for the same response Rg as previously g was, but now surface averaged instead of volume averaged:

 ⌠⌡ surface d2 s  g ·f = ⌠⌡ volume d3 r  gS ·f
Therefore, the line integrals Il,i in the formula given above for track-length estimates, equation 2.2, now reduce to:

 Il,i = ⌠⌡ ri+1 ri dl  δ(ρ) ·g* (l) = g* (ρ = 0) / cos( eρ, Ω)
(2.5)
if the surface is intersected during the test-flight from ri to ri+1 , and Il,i = 0 else. This can easily be seen by writing
 rl,i
 =
ri + l ·

 cos(Ω,eρ)
 cos(Ω, e ρ′ )
 cos(Ω,e ρ")

(2.6)
i.e., ρ = rρ, i + l ·cos( Ω,eρ ) .
Consequently, a statement in subroutine UPSUSR(WEIGHT,IND) for this surface averaged tally must read:
here IS is the index of the surface which is being crossed, ITALS is the labelling number of the surface tally. Furthermore,
 WC = kω*i / [ |v| ·cos( Ω, eρ ) ] ,
g is evaluated at the strike point rS and, as above, kω*i is the weight of the history k after event i and before event i+1. Note:
 WEIGHT = k w*i = k ^ wi+1
in our terminology.
The flag IND has the value 1, if the particle is incident onto the surface, and IND = 2 if the particle is emitted from the surface.
Some care is needed because the local surface normal unit vector eρ at the strike point rS is not known in subroutine UPSUSR unless the surface input variable ILIIN (see input block 3) is equal to 1 or equal to 3 for the surface S. In all other cases, this vector (defined by its Cartesian components CRTX,CRTY,CRTZ in EIRENE) has to be computed in subroutine UPSUSR at each entry, if it is needed for updating the surface averaged tally.
Subroutine UPSUSR is called whenever the default surface tallies are updated. For one sided surface tallies (ILIIN=-2 or ILIIN=-4), therefore, UPSUSR is also only called for one sided tallies only.
In addition to these calls, UPSUSR is also called from surface source routines. For test particles emitted from source surfaces (or generated from incident bulk ions after reflection or sputtering) the call is UPSUSR(WEIGHT,2). For bulk ions (ITYP=4) incident onto a source surface the call is UPSUSR(-WEIGHT,1).
Therefore, the tallies ADDS include the direct source contribution, if the source is a surface source, whereas the default surface tallies don't.

3.3  The user surface reflection model REFUSR

General remarks:

This subroutine is called in the initialization phase of an EIRENE run, at entry RF0USR, from subroutine REFLEC, and at entry SP0USR, from subroutine SPUTER. At these entries the user supplied reflection models and sputtering models (if any), respectively, may be initialized.
Later, during the Monte Carlo sampling phase, it is called via entry RF1USR, whenever a test flight intersects a surface NLLI, for which the fast particle reflection model flag ILREF(NLLI) has been set equal to ILREF(NLLI)=3.
Furthermore, if for a certain surface NLLI a user specified "sputter model" is activated by the flag ILSPT(NLLI)=3, then REFUSR is called at the entry SP1USR, whenever this surface is intersected by a test flight.

Format of subroutine

      SUBROUTINE REFUSR
C USER SUPPLIED SURFACE INTERACTION ROUTINE
C INPUT : CO-ORDINATES OF INCIDENT PARTICLE
C OUTPUT: CO-ORDINATES OF EMITTED PARTICLE
* CALL PARMMOD
:
:
C  INITIALIZE USER SUPPLIED REFLECTION MODEL
ENTRY RF0USR
:
C  DEFINE SPECIES DEPENDENT RECYCLING COEFFICIENT HERE
C     RECYCT(ISPZ,MSURF)=.....
:
RETURN
ENTRY RF1USR (XMW,XCW,XMP,XCP,IGASF,IGAST,ZCOS,ZSIN,EXPI,
.              RPROB,E0TERM,*,*,*,*)
C  RETURN 1: EIRENE STANDARD ANGULAR DISTRIBUTION (DEP. ON EXPI")
C  RETURN 2: THERMAL MOLECULE MODEL (DEP. ON IGAST", E0TERM")
C  RETURN 3: THERMAL ATOM MODEL (DEP. ON IGAST", E0TERM")
C  RETURN 4: ABSORB PARTICLE AT THIS SURFACE
:
:
RETURN
C
C  INITIALIZE USER SUPPLIED SPUTTERING MODEL
ENTRY SP0USR
:
C  DEFINE SPECIES DEPENDENT PHYSICAL SPUTTERING COEFFICIENT HERE
C     RECYCS(ISPZ,MSURF)=.....
C  AND ALSO CHEMICAL SPUTTERING COEFFICIENT
C     RECYCC(ISPZ,MSURF)=.....
:
ENTRY SP1USR
:
:
RETURN
END



3.4  The user source sampling routine SAMUSR

General remarks:

This subroutine is called from the point source sampling routine SAMPNT, the surface source sampling routine SAMSRF, or from the volume sampling routine SAMVOL, if the flag SORLIM(N, ISTRA) (input block 7) for the spatial distribution of birth points on sub-stratum N for stratum ISTRA has a negative value.
This subroutine returns to EIRENE the (cartesian) coordiantes of the birth point
X0,Y0,Z0
cell index information
IRUSR,IPUSR,ITUSR,IAUSR,IBUSR
as well as local plasma background data at this place of birth:
TIWL,TEWL,DIWL,VXWL,VYWL,VZWL,EFWL,SHWL,WEISPZ
The input parameters (EIRENE input block 7 for primary sources)
can be used in this problem specific routine.
Dimensions, Precision:

Format of subroutine

      SUBROUTINE SAMUSR (ISR,X0,Y0,Z0,
.                   IRUSR,IPUSR,ITUSR,IAUSR,IBUSR,
.                   TIWL,TEWL,DIWL,VXWL,VYWL,VZWL,EFWL,SHWL,
.                   WEISPZ)
USE PRECISION
USE PARMMOD
:
:  ! MORE MODULES, IF NEEDED
:
IMPLICIT NONE
INTEGER, INTENT(IN) :: ISR, ISTR
REAL(DP), INTENT(OUT) ::
.          X0,Y0,Z0
INTEGER,  INTENT(OUT) :: IRUSR,IPUSR,ITUSR,IAUSR,IBUSR
!   NPLS: NUMBER OF BULK ION SPECIES IN EIRENE INPUT BLOCK 5
!   NSPZ: TOTAL NUMBER OF PARTICLE SPECIES:
!         NSPZ=NPHOT+NATM+NMOL+NION+NPLS,
!         SEE INPUT BLOCKS 4 AND 5.
REAL(DP), INTENT(OUT) ::
.          TEWL,SHWL,
.          TIWL(NPLS),DIWL(NPLS),
.          VXWL(NPLS),VYWL(NPLS),VZWL(NPLS),
.          EFWL(NPLS),WEISPZ(NSPZ)

!  INITIALIZE SAMPLING FOR SUBSTRATUM  ISR OF STRATUM ISTR
!  ISR:  NUMBER OF SUBSTRATUM (SEE INPUT BLOCK 7)
!  ISTR: NUMBER OF STRATUM

RETURN

ENTRY SM1USR (ISR,X0,Y0,Z0,
.              IRUSR,IPUSR,ITUSR,IAUSR,IBUSR,
.              TIWL,TEWL,DIWL,VXWL,VYWL,VZWL,EFWL,SHWL,WEISPZ)
!  FIND BIRTH POINT COORDINATES. SUBSTRATUM ISR IF STRATUM ISTRA
!  HAS ALREADY BEEN IDENTIFIED. ISTRA (IF NEEDED) IS IN MODULE COMPRT.
!  INITIALIZE THOSE VARIABLES WHICH ARE NOT SET BELOW:
X0  =0._DP
Y0  =0._DP
Z0  =0._DP

IRUSR =0
IPUSR =0
ITUSR =1
IAUSR =0
IBUSR =1

TEWL=0._DP
SHWL=0._DP
TIWL=0._DP
DIWL=0._DP
VXWL=0._DP
VYWL=0._DP
VZWL=0._DP
EFWL=0._DP
WEISPZ=0._DP

!  HERE COMES THE PROBLEM SPECIFIC DEFINITION OF BIRTH POINT SAMPLING

RETURN

END


Note: The parameters EFWL, SHWL have been introduced in 2005.
Note: Correction in May 2006: In previous versions of this manual the parameters TEWL and TIWL have been interchanged with respect to the calling statements in the EIRENE code.
(Thanks to Jose Guasp, (CIEMAT, Spain) for pointing this out.)

to be written

to be written

3.6  The user routines for profiles PROUSR

General remarks:

This subroutine is called from the subroutine PLASMA if the input flag INDPRO (input block 5) has the value INDPRO(IPRO)=5 for a particular background tally IPRO. Only the data for cells from the standard mesh (ICELL=1,..,NSURF) are used. Default "vacuum data" are set for all cells in the "additional cell region" ICELL=NSURF+1,...,NSBOX.
PROUSR is called, furthermore, from subroutine INPUT to provide cell volumes for all cells (ICELL=1,...,NSBOX), if the input flag INDPRO(12) has the value INDPRO(12)=5

Format of subroutine

      SUBROUTINE PROUSR (RHO,INDEX,P0,P1,P2,P3,P4,P5,PVAC,NDAT)
* CALL PARMMOD
:
:
IF (INDEX.EQ.0) THEN
C  DATA FOR TE-PROFILE, NDAT=NSURF
DO ICELL=1,NDAT
RHO(ICELL)=......
ENDDO
ELSEIF (INDEX.EQ.....) THEN
:
:
ELSEIF (INDEX.EQ.5+5*NPLS) THEN
C  DATA FOR VOL, NDAT=NSBOX
:
:
ELSE
WRITE (6,*) 'INVALID INDEX IN PROUSR '
CALL EXIT
ENDIF
RETURN
END


Depending upon the value of INDEX, one has to specify a background tally on the array RHO, RHO(J),J=1,NDAT. These are:
• INDEX = 0 :
Electron temperature (eV)
• INDEX = 1 :
Ion temperature (eV), NPLSI calls, one for each ion species. I.e.:
DO IPLS=1,NPLSI
CALL PROUSR(HELP,1+0*NPLSI,TI0(IPLS),...,NSURF)
CALL RESETP(TIIN,HELP,IPLS,1,NPLS,NSURF)
ENDDO


• INDEX = 1+ NPLS :
Ion density (cm−3), NPLSI calls, one for each ion species
• INDEX = 1+2*NPLS :
Ion drift velocity, x direction (cm/s), NPLSI calls, one for each ion species
• INDEX = 1+3*NPLS :
Ion drift velocity, y direction (cm/s), NPLSI calls, one for each ion species
• INDEX = 1+4*NPLS :
Ion drift velocity, z direction (cm/s), NPLSI calls, one for each ion species
• INDEX = 1+5*NPLS :
Magnetic field, x component, (AU)
• INDEX = 2+5*NPLS :
Magnetic field, y component, (AU)
• INDEX = 3+5*NPLS :
Magnetic field, z component, (AU)
• INDEX = 4+5*NPLS :
Magnetic field, magnitude, (T)
• INDEX = 5+5*NPLS :
Cell volume, (cm3)
• INDEX = 6+5*NPLS :

3.7  User supplied post-processed tally routine TALUSR

to be written
      SUBROUTINE TALUSR(ICOUNT,VECTOR,TALTOT,TALAV,
.                  TXTTL,TXTSP,TXTUN,ILAST,*)



3.8  User supplied "general geometry block"

The following routines have to be provided for the general geometry options (LEVGEO=10 option, NLGEN=TRUE), in which no specific geometry data are available from EIRENE input, and all the geometrical parameters for particle tracing and scoring of tallies are transferred from outside.
• INIUSR (initialize user specified geometry block), see 3.8.1
• LEAUSR(x) (return cell number, for any given position x), see 3.8.2
• TIMUSR (flight time to next cell boundary, and next cell number or surface number), see 3.8.3
• VOLUSR (volume of each grid cell), see 3.8.4
• NORUSR (outer surface normal, for any given position on a surface), see 3.8.5
In case NLGEN, only a reduced set of the standard EIRENE options is available, and the five user routines mentioned above and described below may have to be supplemented by some further options from the problem specific segment USER.F .
Input block 5
only the INDPRO = 3, 5 and 6 options are available. INDPRO = 3 provides constant profiles, same value in each cell. Hence, in case of non-constant background parameters one has to resort to subroutine PROUSR (INDPRO=5), or to the code segment INFCOP (INDPRO=6) for coupling to another source (code) for the data of the background medium.
Input block 7
only the point source option is available. In case of surface sources or volume sources, the subroutine SAMUSR has to be called, see 3.4.
Input block 11
only the printout options are available, no graphics output options are available. However, subroutine PLTUSR(LOG,J) is called from the 2D geometry plotting routine PLT2D (LOG=.TRUE., J = number of surface to be plotted) and from the 3D geometry plotting routine PLT3D (LOG=.FALSE., J = number of surface to be plotted)

3.8.1  Subroutine INIUSR

This subroutine is called from subroutine INPUT, after reading from the formatted input file and after (optional) calls to interfacing routines INFCOP. Any initialization (including definition of COMMON blocks) for the user geometry package may be done here.

3.8.2  Subroutine LEAUSR

The function LEAUSR is called from functions LEARC1 and LEARCA in case of LEVGEO=10. The call is
NCELL=LEAUSR(X,Y,Z)


Here X,Y,Z are the Cartesian coordinates of a point inside the computational domain in centimeters. LEAUSR must then return the number of the cell to which this point belongs.

3.8.3  Subroutine TIMUSR

The subroutine TIMUSR is called from TIMER (block GEO3D). The call is
CALL TIMUSR (NRCELL,X0,Y0,Z0,VELX,VELY,VELZ,
NJUMP,NEWCEL,TIM,ICOS,IERR)


Input
NRCELL
Actual cell number, for which the intersection to the cell boundary has to be found. If NJUMP=0, i.e., for the first call of a new trajectory (e.g., after a collision), the starting point X0,Y0,Z0 must be in that cell. In later calls for the same trajectory (NJUMP > 0), NRCELL has been automatically updated, i.e., it must not necessarily contain the starting point X0,Y0,Z0.
X0,Y0,Z0
Cartesian coordinates of the starting point of a trajectory.
VELX,VELY,VELZ

unit (speed-) vector pointing in the direction of the flight.
Output
NJUMP
= 0
this is the first call for a particular trajectory.
≠ 0
this is a later call for a particular trajectory. The initial position and speed unit vector are the same as in the previous call. Hence: the geometrical parameters depending only on those need not be evaluated again.
NEWCEL
< 0
trajectory has intersected one of the non-default surface (input block 3a). ABS(NEWCEL) is the number of that surface (corresponding to running index ISTSI in input block 3a).
> 0
number of next cell, in which the flight would continue, if no collision event stops the track already earlier. I.e., cell number of the neighbor cell in the direction of the flight.
TIM
distance (cm) from X0,Y0,Z0 to the nearest intersection of the trajectory with a boundary of cell NRCELL
ICOS
only relevant, if the trajectory has intersected a non-default surface. In this case, ICOS is the sign (+1 or -1) of the cosine of the angle of incidence against the surface normal (this later unit vector may, e.g., be found in TIMUSR by a call to NORUSR).
IERR
error flag. Presently any value different from 0 will lead to an immediate end of the run. See subroutine TIMER in code segment GEO3D.F

3.8.4  Subroutine VOLUSR

The subroutine VOLUSR is called from subroutine VOLUME in case of LEVGEO=10. The call is
CALL VOLUSR(N,A)


Here N is the number of cells in this run, and A is an array of length N, containing the N cell volumes in cm3.

3.8.5  Subroutine NORUSR

The subroutine NORUSR is called from subroutine STDUSR in case of LEVGEO=10. The call is
CALL NORUSR(M,X,Y,Z,CX,CY,CZ,SCOS)


Here X,Y,Z are the Cartesian coordinates of a point, located on non-default standard surface no. M. The routine must return the surface normal unit vector CX,CY,CZ.

Chapter 4 Routines for interfacing with other codes: EIRCOP

General remarks

The names of all routines in the interfacing block end with ...COP.

4.1  Routine for interfacing INFCOP

To write an interfacing routine INFCOP in order to couple EIRENE to another code is already a quite formidable task and it is highly recommended that, in addition to the information given below, the user should ask for a few sample routines INFCOP, which are available from the authors.
Data are transported from subroutine INFCOP into EIRENE via the EIRENE work array RWK (Common CSPEI). They are read onto the EIRENE arrays for input tallies (5.1) by calls of subroutine PROFR in the initialization phase, if the flag INDPRO has the value 6 or 7 for a particular input tally (see section 2.5. The following addresses are foreseen on RWK for this data transferring procedure:
Plasma data (INDPRO = 6 option)
TEIN  :  (RWK ((0         ) * NRAD + J), J=1,NSBOX)
TIIN  :  (RWK ((1 + 0*NPLS) * NRAD + J), J=1,NSBOX,I=1,NPLS)
DIIN  :  (RWK ((1 + 1*NPLS) * NRAD + J), J=1,NSBOX,I=1,NPLS)
VXIN  :  (RWK ((1 + 2 NPLS) * NRAD + J), J=1,NSBOX,I=1,NPLS)
VYIN  :  (RWK ((1 + 3*NPLS) * NRAD + J), J=1,NSBOX,I=1,NPLS)
VZIN  :  (RWK ((1 + 4*NPLS) * NRAD + J), J=1,NSBOX,I=1,NPLS)
BXIN  :  (RWK ((1 + 5*NPLS) * NRAD + J), J=1,NSBOX)
BYIN  :  (RWK ((2 + 5*NPLS) * NRAD + J), J=1,NSBOX)
BZIN  :  (RWK ((3 + 5*NPLS) * NRAD + J), J=1,NSBOX)
BFIN  :  (RWK ((4 + 5*NPLS) * NRAD + J), J=1,NSBOX)
VLIN  :  (RWK ((5 + 5*NPLS) * NRAD + J), J=1,NSBOX)


Geometrical data (INDGRD = 6 option)
LEVGEO = 1 or LEVGEO = 2:
RSURF :  (RWK ((6+5*NPLS+NAIN) * NRAD +                 J), J=1,N1ST)
EP    :  (RWK ((6+5*NPLS+NAIN) * NRAD +   N1ST +        J), J=1,N1ST)
EL    :  (RWK ((6+5*NPLS+NAIN) * NRAD + 2*N1ST +        J), J=1,N1ST)
TR    :  (RWK ((6+5*NPLS+NAIN) * NRAD + 3*N1ST +        J), J=1,N1ST)
PSURF :  (RWK ((6+5*NPLS+NAIN) * NRAD + 4*N1ST +        J), J=1,N2ND)
TSURF :  (RWK ((6+5*NPLS+NAIN) * NRAD + 4*N1ST + N2ND + J), J=1,N3RD)
LEVGEO = 3
PGINTF:  (RWK ((6+5*NPLS+NAIN) * NRAD + J), J=1,NPMAX)
TSURF :  (RWK ((6+5*NPLS+NAIN) * NRAD + NPMAX + J), J=1,N3RD)
LEVGEO = 4
TRINTF:  (RWK ((6+5*NPLS+NAIN) * NRAD + J), J=1,NTMAX)
TSURF :  (RWK ((6+5*NPLS+NAIN) * NRAD + NTMAX + J), J=1,N3RD)
LEVGEO = 5
THINTF:  (RWK ((6+5*NPLS+NAIN) * NRAD + J), J=1,NHMAX)

`

4.1.1  entry IF0COP

Geometry data (INDGRD = 6 option)

In the LEVGEO=3 option PGINTF is an array of length NPMAX, which contains all relevant information to generate a 2D mesh of polygons in the x-y plane. It is equivalenced to the common block CPOLYG via the statement:
EQUIVALENCE (PGINTF(1),XPLG(1,1)),...
In the LEVGEO=4 option TRINTF is an array of length NTMAX, which contains all relevant information to generate a mesh of triangles in the x-y plane. It is equivalenced to the common block CTRIA via the statement:
EQUIVALENCE (TRINTF(1),XT(1)),...
In the LEVGEO=5 option THINTF is an array of length NHMAX, which contains all relevant information to generate a mesh of tetrahedrons in the 3d computational domain. It is equivalenced to the common block CTETRA via the statement:
EQUIVALENCE (THINTF(1),XT(1)),...

4.1.2  entry IF1COP

Plasma (background) data (INDPRO = 6 option)

to be written

4.1.3  entry IF2COP(ISTRA)

Primary source data (INDSRC = 6 option)

In case of an EIRENE run, in which input information is obtained from an external code (e.g.: B2, EMC3, DIVIMP, ...), the stratification of the primary source should be such that the first NTARGI (see 2.14) strata are determined by the plasma fluxes onto surfaces ("recycling surface sources"). In this case the source can be defined automatically from the interfacing routine, using the data specified in input block 14. The remaining primary sources NTARGI+1, ..., NSTRAI (e.g. gas puff, volume recombination sources, etc...) still have to defined in input block 7. The input flag INDSRC described below is irrelevant for those strata. Whether input block 7 or input block 14 is used for a particular surface source ISTRA ≤ NTARGI is controlled by the flag INDSRC.
if INDSRC(ISTRA) < 0 then interfacing routine IF2COP is not called. Input data are used from block 7, see 2.7.

if INDSRC(ISTRA) = 0-5, then the input flags described in section 2.7 are used to define the spatial distribution of a recycling source on a grid boundary (function STEP(ISTEP), see section 2.7.1) as well as the distribution in velocity space. The step function itself and the total source strength FLUX(ISTRA), however, are defined in IF2COP, using data from input block 14, see 2.14, for target recycling source ITARG. In this case the sampling range from the spatial surface distribution and the species distribution specified in input block 7 must be equal to or a subrange of the the step-function for target recycling source ITARG defined from the data in input block 14.

if INDSRC(ISTRA) = 6, then the input flags described in section 2.14 are used to define the entire surface recycling source automatically. The corresponding data in input block 7 for this stratum are not used and may be omitted. Stratum ISTRA corresponds to the target recycling source ITARG from block 14.

4.1.4  entry IF3COP(ISTRAA,ISTRAE)

Return data at the end of an EIRENE stratum

At this entry data are transferred back from EIRENE to the external code. This entry is called from the "strata-loop" (subr. MCARLO), after all trajectories for each particular stratum have been sampled and after all volume and surface tallies have been scaled and processed to their final form.
The call to IF3COP is controlled by the flag NMODE (input block 1).
At IF3COP data are expected and prepared for transfer to the external code for strata all ISTRA in the range ISTRA=ISTRAA,ISTRAE.

4.1.5  entry IF4COP

post-processing after one complete cycle: overall balances, etc.

4.2  Routines for cycling of EIRENE with external codes: EIRSRT

In this subroutine the "cycling" between EIRENE and external codes is controlled. There are various options, such as "full time dependent (explicit)", "quasi-stationary (explicit)", and "quasi-stationary (implicit)".

4.3  Routines for special tallies needed for code coupling: UPTCOP

In older versions of EIRENE (before 2002) in particular momentum exchange rates (i.e. friction terms) in the direction parallel to the magnetic field have been programmed here. Meanwhile these have become default tallies, see chapter 5.1. Still some particular coupling tallies, e.g. to render the coupling procedure more implicit, are scored in this routine.

4.4  Statistical noise in Monte Carlo terms for external code, noise-residuals: STATIS_COP

As pointed out in Section 2.9 for all primary (not derived) EIRENE tallies the empirical standard deviation can also be obtained, if requested.
The considerable CPU penalty for doing this was avoided in EIRENE versions 99 and younger, by major code optimization in these parts.
Subroutine STATIS_COP provides these estimates for those tallies which are specific to a particular coupled case. Usually these will be source terms for particle, momentum, and energy balance equations, summed over all donor species.
At entry IF4COP (see section 4.1.5 above) these may be further processed. For example, in case of coupling to the B2 plasma fluid code these noise estimates are integrated into global quantities (dimension: 1/time) and represent the contribution of statistical noise to the overall residuals of a B2 run.
These are printed from IF4COP, together with the overall balances of a coupled run (TRCBAL=.true.)
Hence: if the B2 estimated residuals are of the same order as these "statistical noise residuals", then a further convergence of the combined code can only be achieved by increasing the CPU-time for EIRENE.
Otherwise: The coupled B2-EIRENE run has failed to converge to the solution within the statistical noise.
This is also represented by the convergence measure of "saturated residuals"

Chapter 5 Default EIRENE tallies, and selected Modules

5.1  Tables of EIRENE tallies

The following three tables comprise the EIRENE "tallies". The term "tally" is adopted from neutron transport applications, a more precise terminology would be "response", see section 3.2. Because of the new type of particle introduced during 2002 (photon gas, ITYP=0) the numbering of tallies has changed significantly. Therefore the tables are given below once for the current EIRENE version (2002 and younger), and once for the versions older than 2001.
Each volume averaged response is a spatial function, averaged over a grid cell, hence: piecewise constant in each grid cell ICELL:
Surface averaged responses are also spatial functions, averaged over surfaces, hence: piecewise constant on each surface ISURF or surface segment:
ISURF=1,NLIMI.
and with some spatial resolution possible on standard mesh surfaces
ISURF=NLIM,NLIM+NSTSI
Emitted particle and energy surface averaged fluxes are split by the typ of the incident particle: photons (Phs.), atoms (Ats.), molecules (Mls.), test ions (T.I.) or bulk ions (B.I.). This permits scaling of these fluxes with the factors FPHOT, FATM, FMOL, FION, to eliminate statistical errors from the particle balances (see NLSCL option, input block 1).
If temporal resolution has been requested (input block 13), then NSTSI is increased by one: NSTSIP=NSTSI+1, and the last tally I2=NLIM+NSTSIP corresponds to the "Time-Surface" (the census array) then.

5.1.1  Current status, incl. photon gas tallies (Eirene-02 and younger)

Table 5.1: Input Tallies for Background, Input, Module: COMUSR

 No Name Macroscopic quantity 1.Dim. Units Estim. -1 TEIN Plasma Temperature, Electrons 1 eV / -2 TIIN Plasma Temperature, Bulk Ions NPLS eV / -3 DEIN Plasma Density, Electrons 1 cm−3 / -4 DIIN Plasma Density, Bulk Ions NPLS cm−3 / -5 VXIN Plasma Drift Velocity, x-component, Bulk Ions NPLS cm/sec / -6 VYIN Plasma Drift Velocity, y-component, Bulk Ions NPLS cm/sec / -7 VZIN Plasma Drift Velocity, z-component, Bulk Ions NPLS cm/sec / -8 BXIN Magnetic field unit vector, x-component 1 / / -9 BYIN Magnetic field unit vector, y-component 1 / / -10 BZIN Magnetic field unit vector, z-component 1 / / -11 BFIN Magnetic field strength 1 Tesla / -12 ADIN Additional input tally NAIN see INFCOP / -13 EDRIFT Kinetic energy in drift motion, Bulk Ions NPLS eV / -14 VOL Zone Volume 1 cm3 / -15 WGHT Space and species dependent importance NSPCMC 1 /

Note: DEIN and EDRIFT are "derived" input tallies, internally computed from the ion densities (quasi-neutrality), and the flow field, respectively.
 Table 5.2: Volume Averaged Tallies, Output, Module: CESTIM No Name Macroscopic quantity 1.Dim. Units Estim. Table 5.3: (continued) No Name Macroscopic quantity 1.Dim. Units Estim. 1 PDENA Particle density, Atoms NATM cm−3 T1 2 PDENM Particle density, Molecules NMOL cm−3 T2 3 PDENI Particle density, Test Ions NION cm−3 T3 4 PDENPH Particle density, Photons NPHOT cm−3 T3 5 EDENA Energy density, Atoms NATM eV*cm−3 T4 6 EDENM Energy density, Molecules NMOL eV*cm−3 T5 7 EDENI Energy density, Test Ions NION eV*cm−3 T6 8 EDENPH Energy density, Photons NPHOT eV*cm−3 T6 9 PAEL Particle Source (Electrons) from atom-plasma coll. 1 amp* cm−3 T7 10 PAAT Particle Source (Atoms) from atom-plasma coll. NATM amp*cm−3 T8 11 PAML Particle Source (Molecules) from atom-plasma coll. NMOL amp*cm−3 T9 12 PAIO Particle Source (Test Ions) from atom-plasma coll. NION amp*cm−3 T10 13 PAPHT Particle Source (Photons) from atom-plasma coll. NPHOT amp*cm−3 T10 14 PAPL Particle Source (Bulk Ions) from atom-plasma coll. NPLS amp*cm−3 T11 15 PMEL Particle Source (Electrons) from molecule-plasma coll. 1 amp*cm−3 T12 16 PMAT Particle Source (Atoms) from molecule-plasma coll. NATM amp*cm−3 T13 17 PMML Particle Source (Molecules) from molecule-plasma coll. NMOL amp*cm−3 T14 18 PMIO Particle Source (Test Ions) from molecule-plasma coll. NION amp*cm−3 T15 19 PMPHT Particle Source (Photons) from molecule-plasma coll. NPHOT amp*cm−3 T15 20 PMPL Particle Source (Bulk Ions) from molecule-plasma coll. NPLS amp*cm−3 T16 21 PIEL Particle Source (Electrons) from test ion-plasma coll. 1 amp*cm−3 T17 22 PIAT Particle Source (Atoms) from test ion-plasma coll. NATM amp*cm−3 T18 23 PIML Particle Source (Molecules) from test ion-plasma coll. NMOL amp*cm−3 T19 24 PIIO Particle Source (Test Ions) from test ion-plasma coll. NION amp*cm−3 T20 25 PIPHT Particle Source (Photons) from test ion-plasma coll. NPHOT amp*cm−3 T20 26 PIPL Particle Source (Bulk Ions) from test ion-plasma coll. NPLS amp*cm−3 T21 27 PPHEL Particle Source (Electrons) from photon-plasma coll. 1 amp*cm−3 T17 28 PPHAT Particle Source (Atoms) from photon-plasma coll. NATM amp*cm−3 T18 29 PPHML Particle Source (Molecules) from photon-plasma coll. NMOL amp*cm−3 T19 30 PPHIO Particle Source (Test Ions) from photon-plasma coll. NION amp*cm−3 T20 31 PPHPHT Particle Source (Photons) from photon-plasma coll. NPHOT amp*cm−3 T20 32 PPHPL Particle Source (Bulk Ions) from photon-plasma coll. NPLS amp*cm−3 T21 33 EAEL Energy Source (Electrons) from atom-plasma coll. 1 watt*cm−3 T22 34 EAAT Energy Source (Atoms) from atom-plasma coll. 1 watt*cm−3 T23 35 EAML Energy Source (Molecules) from atom-plasma coll. 1 watt*cm−3 T24 36 EAIO Energy Source (Test Ions) from atom-plasma coll. 1 watt*cm−3 T25 37 EAPHT Energy Source (Photons) from atom-plasma coll. 1 watt*cm−3 T25 38 EAPL Energy Source (Bulk Ions) from atom-plasma coll. 1 watt*cm−3 T26 39 EMEL Energy Source (Electrons) from molecule-plasma coll. 1 watt*cm−3 T27 40 EMAT Energy Source (Atoms) from molecule-plasma coll. 1 watt*cm−3 T28 41 EMML Energy Source (Molecules) from molecule-plasma coll. 1 watt*cm−3 T29 42 EMIO Energy Source (Test Ions) from molecule-plasma coll. 1 watt*cm−3 T30 43 EMPHT Energy Source (Photons) from molecule-plasma coll. 1 watt*cm−3 T30 44 EMPL Energy Source (Bulk Ions) from molecule-plasma coll. 1 watt*cm−3 T31 45 EIEL Energy Source (Electrons) from test ion-plasma coll. 1 watt*cm−3 T32 46 EIAT Energy Source (Atoms) from test ion-plasma coll. 1 watt*cm−3 T33 47 EIML Energy Source (Molecules) from test ion-plasma coll. 1 watt*cm−3 T34 48 EIIO Energy Source (Test Ions) from test ion-plasma coll. 1 watt*cm−3 T35 49 EIPHT Energy Source (Photons) from test ion-plasma coll. 1 watt*cm−3 T35 50 EIPL Energy Source (Bulk Ions) from test ion-plasma coll. 1 watt*cm−3 T36 51 EPHEL Energy Source (Electrons) from photon-plasma coll. 1 watt*cm−3 T33 52 EPHAT Energy Source (Atoms) from photon-plasma coll. 1 watt*cm−3 T33 53 EPHML Energy Source (Molecules) from photon-plasma coll. 1 watt*cm−3 T34 54 EPHIO Energy Source (Test Ions) from photon-plasma coll. 1 watt*cm−3 T35 55 EPHPHT Energy Source (Photons) from photon-plasma coll. 1 watt*cm−3 T35 56 EPHPL Energy Source (Bulk Ions) from photon-plasma coll. 1 watt*cm−3 T36 57 ADDV Additional volume av. Tally, Track-length estimated NADV see UPTUSR TRL 58 COLV Additional volume av. Tally, Collision estimated NCLV see UPCUSR COL 59 SNPV Additional volume av. Tally, Snapshot estimated NSNV see UPSUSR SNP 60 COPV Tallies for coupling to ext. code, see: Subr. UPTCOP NCPV see UPTCOP COP 61 BGKV Volume averaged tallies for BGK-self-collision terms NBGV see BGK T41 62 ALGV Algebraic expression in volume averaged tallies NALV Input ALG 63 PGENA Generation limit, Atoms NATM amp*cm−3 T1 64 PGENM Generation limit, Molecules NMOL amp*cm−3 T1 65 PGENI Generation limit, Test ions NION amp*cm−3 T1 66 PGENPH Generation limit, Photons NPHOT amp*cm−3 T1 67 EGENA dito, Energy flux, Atoms NATM watt*cm−3 T1 68 EGENM dito, Energy flux, Molecules NMOL watt*cm−3 T1 69 EGENI dito, Energy flux, Test ions NION watt*cm−3 T1 70 EGENPH dito, Energy flux, Photons NPHOT watt*cm−3 T1 71 VGENA dito, momentum flux, Atoms NATM as MAPL T1 72 VGENM dito, momentum flux, Molecules NMOL as MMPL T1 73 VGENI dito, momentum flux, Test ions NION as MIPL T1 74 VGENPH dito, momentum flux, Photons NPHOT as MPHPL T1 75 PPAT primary particle sources rate, Atoms NATM amp*cm-3 C1 76 PPML primary particle sources rate, Molecules NMOL amp*cm-3 C1 77 PPIO primary particle sources rate, Test ions NION amp*cm-3 C1 78 PPPHT primary particle sources rate, Photons NPHOT amp*cm-3 C1 79 PPPL primary particle sources rate, Bulk Ions NPLS amp*cm-3 C1 80 EPAT primary energy sources rate, Atoms 1 watt*cm-3 C1 81 EPML primary energy sources rate, Molecules 1 watt*cm-3 C1 82 EPIO primary energy sources rate, Test ions 1 watt*cm-3 C1 83 EPPHT primary energy sources rate, Photons 1 watt*cm-3 C1 84 EPPL primary energy sources rate, Bulk Ions 1 watt*cm-3 C1 85 VXDENA momentum density, x-direction, Atoms NATM g*cm/s*cm-3 T1 86 VXDENM momentum density, x-direction, Molecules NMOL g*cm/s*cm-3 T1 87 VXDENI momentum density, x-direction, Test Ions NION g*cm/s*cm-3 T1 88 VXDENPH momentum density, x-direction, Photons NPHOT g*cm/s*cm-3 T1 89 VYDENA momentum density, y-direction, Atoms NATM g*cm/s*cm-3 T1 90 VYDENM momentum density, y-direction, Molecules NMOL g*cm/s*cm-3 T1 91 VYDENI momentum density, y-direction, Test Ions NION g*cm/s*cm-3 T1 92 VYDENPL momentum density, y-direction, Photons NPHOT g*cm/s*cm-3 T1 93 VZDENA momentum density, z-direction, Atoms NATM g*cm/s*cm-3 T1 94 VZDENM momentum density, z-direction, Molecules NMOL g*cm/s*cm-3 T1 95 VZDENI momentum density, z-direction, Test Ions NION g*cm/s*cm-3 T1 96 VZDENPL momentum density, z-direction, Photons NPHOT g*cm/s*cm-3 T1 97 MAPL parallel momentum source rate, Atoms NATM amp*g*cm/s*cm-3 TPM 98 MMPL parallel momentum source rate, Molecules NMOL amp*g*cm/s*cm-3 TPM 99 MIPL parallel momentum source rate, Test Ions NION amp*g*cm/s*cm-3 TPM 100 MPHPL parallel momentum sources rate, Photons NPHOT amp*g*cm/s*cm-3 TPM
Note: Source means: if the sign is positive, it is a gain for the specified type of particles; if it is negative, it is a loss.

Estimators EIRENE resorts, by default, to tracklength estimators (3.17). Default tallies are updated ("scoring") in routine UPDATE. All default estimators are constant within a cell m, i.e. depend only upon the cell index m but not on the position r in the cell: gt(s) = const(m). Hence they also do not depend on the position s along the track within a cell. UPDATE calls templates UPTUSR, in which any further quantity can be scored by programming the path integral of any function gt(s), see section 3.2.
In some instances still collision estimators (3.15) are employed but we are gradually removing them (and plan keep them in the code only to provide independent checks for code-verification runs.)

TD Tracklength, particle density. gt = 1/v, v=VEL the test particle's velocity. Note that the path-integral of gt along a trajectory in a cell is equal to the time spend by that history in a cell.
TPM Tracklength, parallel (to B-field) momentum source/sink
 Table 5.4: Surface Averaged Tallies, Output, Module: CESTIM No Name Macroscopic quantity 1.Dim. Units Estim. Table 5.5: (continued) No Name Macroscopic quantity 1.Dim. Units Estim. 1 POTAT Particle Flux, incident, Atoms NATM amp T/C 2 PRFAAT Particle Flux, emitted, Ats.⇒Atoms NATM amp T/C 3 PRFMAT Particle Flux, emitted, Mls.⇒Atoms NATM amp T/C 4 PRFIAT Particle Flux, emitted, T.I.⇒Atoms NATM amp T/C 5 PRFPHAT Particle Flux, emitted, Pht.⇒Atoms NATM amp T/C 6 PRFPAT Particle Flux, emitted, B.I.⇒Atoms NATM amp T/C 7 POTML Particle Flux, incident, Molecules NMOL amp T/C 8 PRFAML Particle Flux, emitted, Ats.⇒Molecules NMOL amp T/C 9 PRFMML Particle Flux, emitted, Mls.⇒Molecules NMOL amp T/C 10 PRFIML Particle Flux, emitted, T.I.⇒Molecules NMOL amp T/C 11 PRFPHML Particle Flux, emitted, Pht.⇒Molecules NMOL amp T/C 12 PRFPML Particle Flux, emitted, B.I.⇒Molecules NMOL amp T/C 13 POTIO Particle Flux, incident, Test Ions NION amp T/C 14 PRFAIO Particle Flux, emitted, Ats.⇒Test Ions NION amp T/C 15 PRFMIO Particle Flux, emitted, Mls.⇒Test Ions NION amp T/C 16 PRFIIO Particle Flux, emitted, T.I.⇒Test Ions NION amp T/C 17 PRFPHIO Particle Flux, emitted, Pht.⇒Test Ions NION amp T/C 18 PRFPIO Particle Flux, emitted, B.I.⇒Test Ions NION amp T/C 19 POTPHT Particle Flux, incident, Photons NPHOT amp T/C 20 PRFAPHT Particle Flux, emitted, Ats.⇒Photons NPHOT amp T/C 21 PRFMPHT Particle Flux, emitted, Mls.⇒Photons NPHOT amp T/C 22 PRFIPHT Particle Flux, emitted, T.I.⇒Photons NPHOT amp T/C 23 PRFPHPHT Particle Flux, emitted, Pht.⇒Photons NPHOT amp T/C 24 PRFPPHT Particle Flux, emitted, B.I.⇒Photons NPHOT amp T/C 25 POTPL Particle Flux, incident, Bulk Ions NPLS amp T/C 26 EOTAT Energy Flux, incident, Atoms NATM watt T/C 27 ERFAAT Energy Flux, emitted, Ats.⇒Atoms NATM watt T/C 28 ERFMAT Energy Flux, emitted, Mls.⇒Atoms NATM watt T/C 29 ERFIAT Energy Flux, emitted, T.I.⇒Atoms NATM watt T/C 30 ERFPHAT Energy Flux, emitted, Pht.⇒Atoms NATM watt T/C 31 ERFPAT Energy Flux, emitted, B.I.⇒Atoms NATM watt T/C 32 EOTML Energy Flux, incident, Molecules NMOL watt T/C 33 ERFAML Energy Flux, emitted, Ats.⇒Molecules NMOL watt T/C 34 ERFMML Energy Flux, emitted, Mls.⇒Molecules NMOL watt T/C 35 ERFIML Energy Flux, emitted, T.I.⇒Molecules NMOL watt T/C 36 ERFPHML Energy Flux, emitted, Pht.⇒Molecules NMOL watt T/C 37 ERFPML Energy Flux, emitted, B.I.⇒Molecules NMOL watt T/C 38 EOTIO Energy Flux, incident, Test Ions NION watt T/C 39 ERFAIO Energy Flux, emitted, Ats.⇒Test Ions NION watt T/C 40 ERFMIO Energy Flux, emitted, Mls.⇒Test Ions NION watt T/C 41 ERFIIO Energy Flux, emitted, T.I.⇒Test Ions NION watt T/C 42 ERFPHIO Energy Flux, emitted, Pht.⇒Test Ions NION watt T/C 43 ERFPIO Energy Flux, emitted, B.I.⇒Test Ions NION watt T/C 44 EOTPHT Energy Flux, incident, Photons NPHOT watt T/C 45 ERFAPHT Energy Flux, emitted, Ats.⇒Photons NPHOT watt T/C 46 ERFMPHT Energy Flux, emitted, Mls.⇒Photons NPHOT watt T/C 47 ERFIPHT Energy Flux, emitted, T.I.⇒Photons NPHOT watt T/C 48 ERFPHPHT Energy Flux, emitted, Pht.⇒Photons NPHOT watt T/C 49 ERFPPHT Energy Flux, emitted, B.I.⇒Photons NPHOT watt T/C 50 EOTPL Energy Flux, incident, Bulk Ions NPLS watt T/C 51 SPTAT Sputtered Flux, by incident Atoms NATM amp T/C 52 SPTML Sputtered Flux, by incident Molecules NMOL amp T/C 53 SPTIO Sputtered Flux, by incident Test Ions NION amp T/C 54 SPTPHT Sputtered Flux, by incident Photons NPHOT amp T/C 55 SPTPL Sputtered Flux, by incident Bulk Ions NPLS amp T/C 56 SPTTOT Sputtered Flux, total 1 amp T/C 57 ADDS Additional Surface Tally NADS Input T/C 58 ALGS Algebraic expression in surface averaged tallies NALS Input 59 SPUMP Pumped flux, by species NSPZ amp
Note: The tallies listed here are two dimensional arrays. The 2nd index I2 is the number of the surface or the surface segment. For I2=1,...,NLIMI the tallies correspond to the "Additional Surfaces", (input block 3B).
For I2=NLIM+1,NLIM+NSTSI*NGITT the data correspond to the "Non-default Standard Surfaces" (input block 3A). On each such surface of block 3A, there is a spatial resolution with up to NGITT surface segments, depending upon the standard grid dimensionality.

5.1.2  old version, w/o photon gas tallies (Eirene-01 and older)

Table 5.6: Input Tallies for Background, Input, Common: COMUSR

 No Name Macroscopic quantity 1.Dim. Units Estim. -1 TEIN Plasma Temperature, Electrons 1 eV / -2 TIIN Plasma Temperature, Bulk Ions NPLS eV / -3 DEIN Plasma Density, Electrons 1 cm−3 / -4 DIIN Plasma Density, Bulk Ions NPLS cm−3 / -5 VXIN Plasma Drift Velocity, x-component, Bulk Ions NPLS cm/sec / -6 VYIN Plasma Drift Velocity, y-component, Bulk Ions NPLS cm/sec / -7 VZIN Plasma Drift Velocity, z-component, Bulk Ions NPLS cm/sec / -8 BXIN Magnetic field unit vector, x-component 1 / / -9 BYIN Magnetic field unit vector, y-component 1 / / -10 BZIN Magnetic field unit vector, z-component 1 / / -11 BFIN Magnetic field strength 1 Tesla / -12 ADIN Additional input tally NAIN see INFCOP / -13 EDRIFT Kinetic energy in drift motion, Bulk Ions NPLS eV / -14 VOL Zone Volume 1 cm3 / -15 WGHT Space and species dependent importance NSPCMC 1 /
Table 5.7: Volume Averaged Tallies, Output, Common: CESTIM
 No Name Macroscopic quantity 1.Dim. Units Estim. 1 PDENA Particle density, Atoms NATM cm−3 T1 2 PDENM Particle density, Molecules NMOL cm−3 T2 3 PDENI Particle density Test Ions NION cm−3 T3 4 EDENA Energy density, Atoms NATM eV*cm−3 T4 5 EDENM Energy density, Molecules NMOL eV*cm−3 T5 6 EDENI Energy density, Test Ions NION eV*cm−3 T6 7 PAEL Particle Source (Electrons) from atom-plasma coll. 1 amp* cm−3 T7 8 PAAT Particle Source (Atoms) from atom-plasma coll. NATM amp*cm−3 T8 9 PAML Particle Source (Molecules) from atom-plasma coll. NMOL amp*cm−3 T9 10 PAIO Particle Source (Test Ions) from atom-plasma coll. NION amp*cm−3 T10 11 PAPL Particle Source (Bulk Ions) from atom-plasma coll. NPLS amp*cm−3 T11 12 PMEL Particle Source (Electrons) from molecule-plasma coll. 1 amp*cm−3 T12 13 PMAT Particle Source (Atoms) from molecule-plasma coll. NATM amp*cm−3 T13 14 PMML Particle Source (Molecules) from molecule-plasma coll. NMOL amp*cm−3 T14 15 PMIO Particle Source (Test Ions) from molecule-plasma coll. NION amp*cm3 T15 16 PMPL Particle Source (Bulk Ions) from molecule-plasma coll. NPLS amp*cm−3 T16 17 PIEL Particle Source (Electrons) from test ion-plasma coll. 1 amp*cm−3 T17 18 PIAT Particle Source (Atoms) from test ion-plasma coll. NATM amp*cm−3 T18 19 PIML Particle Source (Molecules) from test ion-plasma coll. NMOL amp*cm−3 T19 20 PIIO Particle Source (Test Ions) from test ion-plasma coll. NION amp*cm−3 T20 21 PIPL Particle Source (Bulk Ions) from test ion-plasma coll. NPLS amp*cm−3 T21 22 EAEL Energy Source (Electrons) from atom-plasma coll. 1 watt*cm−3 T22 23 EAAT Energy Source (Atoms) from atom-plasma coll. 1 watt*cm−3 T23 24 EAML Energy Source (Molecules) from atom-plasma coll. 1 watt*cm−3 T24 25 EAIO Energy Source (Test Ions) from atom-plasma coll. 1 watt*cm−3 T25 26 EAPL Energy Source (Bulk Ions) from atom-plasma coll. 1 watt*cm−3 T26 27 EMEL Energy Source (Electrons) from molecule-plasma coll. 1 watt*cm−3 T27 28 EMAT Energy Source (Atoms) from molecule-plasma coll. 1 watt*cm−3 T28 29 EMML Energy Source (Molecules) from molecule-plasma coll. 1 watt*cm−3 T29 30 EMIO Energy Source (Test Ions) from molecule-plasma coll. 1 watt*cm−3 T30 31 EMPL Energy Source (Bulk Ions) from molecule-plasma coll. 1 watt*cm−3 T31 32 EIEL Energy Source (Electrons) from test ion-plasma coll. 1 watt*cm−3 T32 33 EIAT Energy Source (Atoms) from test ion-plasma coll. 1 watt*cm−3 T33 34 EIML Energy Source (Molecules) from test ion-plasma coll. 1 watt*cm−3 T34 35 EIIO Energy Source (Test Ions) from test ion-plasma coll. 1 watt*cm−3 T35 36 EIPL Energy Source (Bulk Ions) from test ion-plasma coll. 1 watt*cm−3 T36 37 ADDV Additional volume av. Tally, Track-length estimated NADV see UPTUSR TRL 38 COLV Additional volume av. Tally, Collision estimated NCLV see UPCUSR COL 39 SNPV Additional volume av. Tally, Snapshot estimated NSNV see UPSUSR SNP 40 COPV Tallies for coupling to ext. code, see: Subr. UPTCOP NCPV see UPTCOP COP 41 BGKV Volume averaged tallies for BGK-self-collision terms NBGV see BGK T41 42 ALGV Algebraic expression in volume averaged tallies NALV Input ALG 43 PGENA Generation limit, Atoms NATM amp*cm−3 T1 44 PGENM Generation limit, Molecules NMOL amp*cm−3 T1 45 PGENI Generation limit, Test ions NION amp*cm−3 T1 46 EGENA dito, Energy flux, Atoms NATM watt*cm−3 T1 47 EGENM dito, Energy flux, Molecules NMOL watt*cm−3 T1 48 EGENI dito, Energy flux, Test ions NION watt*cm−3 T1 49 VGENA dito, momentum flux, Atoms NATM as COP T1 50 VGENM dito, momentum flux, Molecules NMOL as COP T1 51 VGENI dito, Momentum flux, Test ions NION as COP T1

Note: Source means: if the sign is positive, it is a gain for the specified type of particles; if it is negative, it is a loss.
Table 5.8: Surface Averaged Tallies, Output, Common: CESTIM

 No Name Macroscopic quantity 1.Dim. Units Estim. 1 POTAT Particle Flux, incident, Atoms NATM amp T/C 2 PRFAAT Particle Flux, emitted, Ats.⇒Atoms NATM amp T/C 3 PRFMAT Particle Flux, emitted, Mls.⇒Atoms NATM amp T/C 4 PRFIAT Particle Flux, emitted, T.I.⇒Atoms NATM amp T/C 5 PRFPAT Particle Flux, emitted, B.I.⇒Atoms NATM amp T/C 6 POTML Particle Flux, incident, Molecules NMOL amp T/C 7 PRFAML Particle Flux, emitted, Ats.⇒Molecules NMOL amp T/C 8 PRFMML Particle Flux, emitted, Mls.⇒Molecules NMOL amp T/C 9 PRFIML Particle Flux, emitted, T.I.⇒Molecules NMOL amp T/C 10 PRFPML Particle Flux, emitted, B.I.⇒Molecules NMOL amp T/C 11 POTIO Particle Flux, incident, Test Ions NION amp T/C 12 PRFAIO Particle Flux, emitted, Ats.⇒Test Ions NION amp T/C 13 PRFMIO Particle Flux, emitted, Mls.⇒Test Ions NION amp T/C 14 PRFIIO Particle Flux, emitted, T.I.⇒Test Ions NION amp T/C 15 PRFPIO Particle Flux, emitted, B.I.⇒Test Ions NION amp T/C 16 POTPL Particle Flux, incident, Bulk Ions NPLS amp T/C 17 EOTAT Energy Flux, incident, Atoms NATM watt T/C 18 ERFAAT Energy Flux, emitted, Ats.⇒Atoms NATM watt T/C 19 ERFMAT Energy Flux, emitted, Mls.⇒Atoms NATM watt T/C 20 ERFIAT Energy Flux, emitted, T.I.⇒Atoms NATM watt T/C 21 ERFPAT Energy Flux, emitted, B.I.⇒Atoms NATM watt T/C 22 EOTML Energy Flux, incident, Molecules NMOL watt T/C 23 ERFAML Energy Flux, emitted, Ats.⇒Molecules NMOL watt T/C 24 ERFMML Energy Flux, emitted, Mls.⇒Molecules NMOL watt T/C 25 ERFIML Energy Flux, emitted, T.I.⇒Molecules NMOL watt T/C 26 ERFPML Energy Flux, emitted, B.I.⇒Molecules NMOL watt T/C 27 EOTIO Energy Flux, incident, Test Ions NION watt T/C 28 ERFAIO Energy Flux, emitted, Ats.⇒Test Ions NION watt T/C 29 ERFMIO Energy Flux, emitted, Mls.⇒Test Ions NION watt T/C 30 ERFIIO Energy Flux, emitted, T.I.⇒Test Ions NION watt T/C 31 ERFPIO Energy Flux, emitted, B.I.⇒Test Ions NION watt T/C 32 EOTPL Energy Flux, incident, Bulk Ions NPLS watt T/C 33 SPTAT Sputtered Flux, by incident Atoms NATM amp T/C 34 SPTML Sputtered Flux, by incident Molecules NMOL amp T/C 35 SPTIO Sputtered Flux, by incident Test Ions NION amp T/C 36 SPTPL Sputtered Flux, by incident Bulk Ions NPLS amp T/C 37 SPTTOT Sputtered Flux, total 1 amp T/C 38 ADDS Additional Surface Tally NADS Input T/C 39 ALGS Algebraic expression in surface averaged tallies NALS Input 40 SPUMP Pumped flux, by species NSPZ amp

Note: The tallies listed here are two dimensional arrays. The 2nd index I2 is the number of the surface or the surface segment. For I2=1,...,NLIMI the tallies correspond to the "Additional Surfaces", (input block 3B).
For I2=NLIM+1,NLIM+NSTSI*NGITT the data correspond to the "Non-default Standard Surfaces" (input block 3A). On each such surface, there is a spatial resolution with up to NGITT surface segments.

Bibliography

[1]
B. Braams. Computational Studies in Tokamak Equilibrium and Transport. PhD thesis, Rijksuniversiteit Utrecht, June 1986.
[2]
D. Reiter, H. Kever, G.H. Wolf, et al. Helium removal from tokamks. Plasma Phys. and Contr. Fus., 33:1579, 1991.
[3]
D. Reiter. Progress in 2-dimensional plasma edge modelling. J. Nucl. Mat., 196-198:241, 1992.
[4]
D. Reiter. The EIRENE code, Version Jan. 92, User manual, March 1992.
[5]
D. Reiter, P. Börner, B. Küppers, M. Baelmans, and G. Maddison. Final report on KFA-NET contract 428/90-8/FU-D (1991).
[6]
G.P. Maddison, E.S. Hotston, D. Reiter, et al. Towards fully authentic modelling of iter divertor plasma's. In Proc. 18th Eur. Conf. on Contr. Fus. and Plasma Phys., volume 15C, page 197, Berlin, 1991.
[7]
J. Spanier and E.M. Gelbard. Monte Carlo Principles and neutron transport problems. Addison Wesley Publication Company, 1969.
[8]
D. Reiter, Chr. May, et al. J. Nucl. Mat., 220:987, 1994.
[9]
G.P. Maddison and D. Reiter. Recycling source terms for edge plasma fluid models and impact on convergence behaviour in the braams b2 code. KFA-Jülich Report Jül- 2872, Forschungszentrum Jülich, March 1994.
[10]
D.B. Heifetz, D. Post, M. Petravic, et al. A monte carlo model of neutral particle transport in diverted plasmas. Princeton report PPPL 1843, PPPL, November 1981. J.Comput.Phys. 46, 309 (1982).
[11]
E. Cupini, A. de Matteis, and R.Simonini. EUR XII-324/9, April 1983.
[12]
L. Devroye. Non-Uniform Random Variate Generation. Springer Verlag, 1986.
[13]
R. Behrisch. Plasma-wall interaction. In Summer School of Tokamak Reactors for Breakeven, Erice, 1976.
[14]
Impurity Control, INTOR Workshop Phase IIA, Part 3.
[15]
W. Eckstein and D.B. Heifetz. Data sets for hydrogen reflection and their use in neutral transport calculations. MPI-Garching Report IPP 9/59, MPI-Garching, August 1986. J.Nucl.Mater. 145-147, p332 (1987).
[16]
G. Bateman. Distribution of neutrals scattered off a wall. PPPL Appl. Phys. Rep. No. 1, PPPL, 1980.
[17]
D. Reiter, P. Bogen, and U. Samm. J. Nucl. Mat., 196-198:1059, 1992.
[18]
D. Reiter, Chr. May, M. Baelmans, et al. J. Nucl. Mat., 241-243:342, 1996.
[19]
Th. Behringer. Einfluß nichtlinearer effekte auf den neutralgastransport in tokamaks. KFA-Jülich Report Jül- 2637, Forschungszentrum Jülich, June 1992. Dissertation.
[20]
W.D. Langer. Nuclear Fusion, 22(6):751, 1982.
[21]
D.L. Book. NRL Plasma Formulary. NRL Publication 0084-4040, Washington, DC 20375, 1987.
[22]
D. Reiser and D. Reiter. Nuclear Fusion, xx:xxx.
[23]
R.K. Janev, W.D. Langer, K. Evans, Jr., et al. Elementary Processes in Hydrogen-Helium Plasmas, volume 4 of Springer Series on Atoms + Plasmas. Springer Verlag, 1987.
[24]
A.B. Ehrhardt and W.D. Langer. Collisional processes of hydrocarbons in hydrogen plasmas. Princeton report PPPL 2477, PPPL, September 1987.
[25]
D. Reiter. Atomic and Plasma-Material Interaction Processes in Controlled Thermonuclear Fusion. Elsevier Science Publishers, 1993.
[26]
Gordeev et al. Pis'ma Zh. Ehksp. Teor. Fiz., (25):223, 1977.
[27]
A. Nicolai and D. Reiter. J. Comp. Phys., 55(1):129-153, 1984.
[28]
Nuclear Fusion, Special issue 1984, Data Compendium for Plasma-Surface interactions.
[29]
W. Eckstein, Garcia-Rosales C., J. Roth, et al. Sputtering data. MPI-Garching Report IPP-9/82, MPI-Garching, February 1993.
[30]
J. Roth and C. Gracia-Rosales. Nucl.Fus. 36, 12:1647, 1996.
[31]
J. Roth. J.Nucl.Mater., 266-269:51-57, 1999.

File translated from TEX by TTHgold, version 3.82.
On 28 May 2010, 11:06.