3-NIRVANA-user-guide/3.2-User-interfaces.md

-![logo_short_small](uploads/1e4ee858ca174663b210647640da0976/logo_short_small.png) `4.1`
+![logo](uploads/1e4ee858ca174663b210647640da0976/logo_short_small_circle.png) `4.2` | Udo Ziegler, AIP
 
 CONTENTS:  
-[Overview](3.2-User-interfaces#overview)  
-[Specification of main simulation parameters](3.2-User-interfaces#specification-of-main-simulation-parameters)    
-[Defining initial conditions](3.2-User-interfaces#defining-initial-conditions)  
-[Defining boundary conditions](3.2-User-interfaces#defining-boundary-conditions)  
-[User-defined coefficients for dissipative processes](3.2-User-interfaces#user-defined-coefficients-for-dissipative-processes)  
-[User-defined coefficient for ambipolar diffusion](3.2-User-interfaces#user-defined-coefficient-for-ambipolar-diffusion)  
-[User-defined body force](3.2-User-interfaces#user-defined-body-force)  
-[User-defined cooling and heating function](3.2-User-interfaces#user-defined-cooling-and-heating-function)  
-[User-defined equation of state](3.2-User-interfaces#user-defined-equation-of-state)  
-[User-defined mesh refinement](3.2-User-interfaces#user-defined-mesh-refinement)  
-[User-defined background magnetic field](3.2-User-interfaces#user-defined-background-magnetic-field)  
-[Specification of NCCM parameters](3.2-User-interfaces#specification-of-nccm-parameters)  
-[User-controllable macros](3.2-User-interfaces#user-controllable-macros)
-
-## Overview
-
-User interfaces are certain code modules in order to generate and
-customize a simulation problem. There are three different types of user
-interfaces in NIRVANA:
-
-\(1\) User interfaces which demand to add or modify code. These are
-recognized by C files ending on `User.c`:
-
--   `configUser.c` – definition of IC
+[Overview](#overview)  
+[Specification of main simulation parameters](#specification-of-main-simulation-parameters)  
+[Defining initial conditions](#defining-initial-conditions)  
+[Defining boundary conditions](#defining-boundary-conditions)   
+[User-defined coefficients for dissipative processes](#user-defined-coefficients-for-dissipative-processes)  
+[User-defined coefficient for ambipolar diffusion](#user-defined-coefficient-for-ambipolar-diffusion)  
+[User-defined body force](#user-defined-body-force)  
+[User-defined cooling- and heating function](#user-defined-cooling-and-heating-function)  
+[User-defined equation of state](#user-defined-equation-of-state)  
+[User-defined initial mesh refinement and refinement control](#user-defined-initial-mesh-refinement-and-refinement-control)   
+[User-defined background magnetic field](#user-defined-background-magnetic-field)   
+[User-defined MHD driver](#user-defined-mhddriver)   
+[Specification of NCCM parameters](#specification-of-nccm-parameters)  
+[User-controllable macros](#user-controllable-macros)   
+
+### Overview
+
+User interfaces are code modules in order to generate and customize a
+simulation problem. There are three different types of user interfaces
+in NIRVANA:
+
+**(I)** User interfaces which demand to add or modify code. These are
+recognized by files ending on `User.c`:
+
+-   `configUser.c` -- definition of IC
 
 -   `bcrhoUser.c`,`bcmUser.c`,`bceUser.c`,`bcetUser.c`,
-    `bcbUser.c`,`bcnXUser.c`,`bcCUser.c`,`bctfUser.c` – user-defined BC
+    `bcbUser.c`,`bcnXUser.c`,`bcCUser.c`,`bctfUser.c` -- user-defined BC
 
 -   `viscosityCoeffUser.c`,`conductionCoeffUser.c`,`diffusionCoeffUser.c`
-    – user-defined dissipation coefficients for fluid viscosity, thermal
-    conduction and magnetic diffusion
+    -- user-defined dissipation coefficients for fluid viscosity,
+    thermal conduction and Ohmic diffusion
 
--   `APdiffusionCoeffUser.c` – user-defined coefficient for ambipolar
+-   `APdiffusionCoeffUser.c` -- user-defined coefficient for ambipolar
     diffusion
 
--   `forceUser.c` – user-defined body force
+-   `forceUser.c` -- user-defined body force
 
--   `sourceCoolingUser.c`,`sourceHeatingUser.c` – user-defined
-    cooling/heating function
+-   `sourceCoolingUser.c`,`sourceHeatingUser.c` -- user-defined cooling-
+    and heating function
 
--   `eosTabPressureUser.c`,`eosTabTemperatureUser.c` – user-defined
+-   `eosTabPressureUser.c`,`eosTabTemperatureUser.c` -- user-defined
     pressure and temperature for a tabulated equation of state
 
--   `initDomainUser.c`,`checkDomainUser.c` – user-defined initial- and
-    restricted mesh refinement
+-   `initDomainUser.c`,`checkDomainUser.c` -- user-defined initial mesh
+    refinement and refinement control
 
--   `sourceB0User.c` - user-defined background magnetic field in the 
+-   `sourceB0User.c` -- user-defined background magnetic field in the
     B-splitting formalism
 
-Defining interfaces requires some familiarity with NIRVANA and, in
+-   `driverMHDUser.c` -- user-defined driver for MHD flows
+
+Coding those interfaces requires some familiarity with NIRVANA and, in
 particular, a basic understanding of its data structure. The defintion
-of IC in `configUser.c` is the minimum necessary for a problem setup.
-Sample `...User.c` files can be found in the various testproblems
-subdirectories under top-level `/nirvana/testproblems`.
+of IC in module `configUser.c` is the least necessary for a problem
+setup. Some examples of `...User.c` files can be found in the
+testproblems subdirectories below `/nirvana/testproblems`.
 
-\(2\) User parameter files (suffix `.par`) which demand editing in order
-to configure simulation properties:
+**(II)** User parameter files (ending on `.par`) which demand editing to
+configure simulation properties:
 
--   `nirvana.par` – specification of main simulation parameters
+-   `nirvana.par` -- specification of main simulation parameters
 
--   `NCCM.par` – specification of NCCM parameters
+-   `NCCM.par` -- specification of NCCM parameters
 
-Configuration of the main parameter file `nirvana.par` is the minimum
-necessary for a problem setup.
+Customization of the parameter file `nirvana.par` is the least necessary
+to configure a problem.
 
-\(3\) User header files ending on `User.h`:
+**(III)** User header files ending on `User.h`:
 
--   `nirvanaUser.h` – header file containing macros for
-    user-controllable, secondary simulation parameters
+-   `nirvanaUser.h` -- header file containing user-contollable macros
 
--   `User.h` – freely available header file
+-   `User.h` -- freely available header file
 
-## Specification of main simulation parameters
+### Specification of main simulation parameters
 
-NIRVANA requires the parameter file `nirvana.par` as input. In this file
-the user specifies main parameters which characterize the properties of
-a simulation. Within `nirvana.par` parameters are grouped into blocks
-with category names
+`nirvana.par` is the place to specify important simulation parameters.
+The file is divided into several sections grouping parameters into
+categories:
 
--   **SIMULATION I/O** – I/O-related parameter
+-   **SIMULATION I/O** -- I/O-related parameter
 
--   **GEOMETRY** – selection of the coordinate system
+-   **GEOMETRY** -- selection of the coordinate system
 
--   **DOMAIN SETTINGS** – physical/numerical domain specifications
+-   **DOMAIN SETTINGS** -- physical/numerical domain specifications
 
--   **BOUNDARY CONDITIONS** – selection of boundary conditions
+-   **BOUNDARY CONDITIONS** -- selection of boundary conditions
 
--   **MESH REFINEMENT** – AMR-related parameters
+-   **MESH REFINEMENT** -- AMR-related parameters
 
--   **USER-SPECIFIC PARAMETERS** – freely usuable parameter
+-   **USER-SPECIFIC PARAMETERS** -- freely usuable parameter
 
--   **SOLVER SPECIFICATIONS** – selection of solver options
+-   **SOLVER SPECIFICATIONS** -- selection of solver options
 
--   **PHYSICS SPECIFICATIONS** – selection of physics options
+-   **PHYSICS SPECIFICATIONS** -- selection of physics options
 
-A parameter section starts with its category name preceded by the
-`>`-character. The associated parameters follow in lines where a line
-can contain more than one parameter. Parameters are separated by white
-spaces and it can be numbers, single characters or words. Text following
-the `>`-character is ignored by the parser and serves for comments. For
-example, the GEOMETRY section in `nirvana.par` looks like
+A parameter section starts with its category name preceded by a
+`>`-character. The parameters follow in lines where a line can contain
+more than one parameter. Parameters are separated by white spaces and it
+can be numbers, single characters or words. Text following the
+`>`-character at the end of a line is comment and is ignored by the
+parser.
 
-    >GEOMETRY ----------------------------------------------------------------------
-     CART  0.00e+00 0.00e+00 0.00e+00            >geometry:{CART,CYL,SPH},omega[0-2]
-
-There is just one parameter line here containing 4 parameters: the first
-specifies the type of coordinate system being `CART` for Cartesian and
-the remaining three parameters give the vector components (doubles) of
-the angular velocity vector of a rotating frame of reference with
-respect to the inertial frame. It follows a comment after `>`.
-
-A complete list of parameters for each category is presented next.
-Parameter lines are marked with a 2-digit running number for a clear
-reference (Lines are actually not numbered in `nirvana.par`). Parameter
-explanations follow the pattern ’line\# (parameter names within the line
-from left to right)’ plus a description of each parameter appearing in a
-line. For some parameters their valid range of values are given in round
-brackets after the parameter’s name either in form of a set {} or a
-numeric range. For some parameters typically used values are given.
-Technically, most of the parameters are collected in a global structure
-variable named `_C` defined in the header `nirvana.h`. A certain
-variable `V` in `_C` is then addressed by `_C.V`, e.g., `_C.geometry`
-for the type of coordinate system.
+All parameters are described in detail in the following. On referencing
+purposes parameter lines are marked here with a 2-digit running number.
+Parameter names set in parentheses follow after. Actually, line numbers
+and names are not present in `nirvana.par`. Note that almost all
+parameters in `nirvana.par` are subsumed in a global structure variable
+named `_C` defined in the header file `ctrl.h`. A certain parameter `V`
+in `_C` is then addressed by `_C.V`, e.g., `_C.data_dir`, e.g.,
+`_C.geometry` for the type of coordinate system. For some parameters
+their range of validity is given either in the form of a set {} of
+possible values or a numeric range.
 
 **SIMULATION I/O**:
 
@@ -136,59 +126,52 @@ for the type of coordinate system.
 
     -   `mode` ({NEW,RES,MOD,ANA}): running mode of a simulation. Value
         NEW starts a new simulation. Value RES restarts a simulation
-        from a snapshot with name specified by the parameter `fname`.
-        Value MOD is like RES but, in addition, allows a modification of
-        read-in simulation data by calling the user interface function
-        `modifyConfigUser()` before the simulation continues. Value ANA
-        reads data from a snapshot followed by a call of the user
-        interface function `analysisUser()` and exits thereafter. The
-        ANA mode is intended for a one-time user-specific analysis of
-        the snapshot data with full availability of all code
-        functionalities. In the modes RES, MOD and ANA all parameters
-        following `fname` in `nirvana.par` are ignored.
-
-    -   `fname`: qualified name of the snapshot. The qualified name here
-        is the path relative to the `NIRVANA_RUN` directory as set in
-        the `runNIRVANA` (or `runNIRVANA_MPI.PBS`) script plus the
-        *base* container filename of the snapshot. The base container
-        filename is `NIR#` or `NIRLAST#` with `#`=`_C.mod` but
-        excluding the dot and post-dot container id (cf.
-        \[Output data\](3.3-Output-data\#snapshot-files)).
+        from a snapshot with filename given by the parameter `fname`.
+        Value MOD is like RES but, in addition, allows modification of
+        the imported data via coding the user interface function
+        `modifyConfigUser()`. Value ANA reads data from a snapshot and
+        then calls the user interface function `analysisUser()` before
+        the code exits. This mode is intended for data analysis of a
+        snapshot with all the code functionality available as in a
+        regular run. In mode RES, MOD and ANA all remaining parameters
+        in `nirvana.par` following `fname` are ignored.
+
+    -   `fname`: snapshot filename including the directory path relative
+        to the `NIRVANA_RUN` directory as set in `runNIRVANA` (or
+        `runNIRVANA_PBS`, `runNIRVANA_SLURM`). The base filename has to
+        be specified without container id file extension, i.e., `NIR#`
+        or `NIRLAST` instead of `NIR#.#` or `NIRLAST.#` (cf. [Snapshot
+        files](#snapshot-files)).
 
 -   `02` (`_C.mod_max`, `_C.time_max`)
 
-    -   `_C.mod_max`: maximum number of timestep cycles a simulation
-        should run. If the physical evolution time reaches the value
-        `_C.time_max` before `_C.mod_max` cycles are done the simulation
-        stops ahead of schedule.
+    -   `_C.mod_max`: maximum number of timesteps a simulation should
+        run. If the physical evolution time reaches the value
+        `_C.time_max` before `_C.mod_max` timesteps are done the
+        simulation stops ahead of schedule.
 
     -   `_C.time_max`: maximum physical evolution time a simulation
-        should run. If `_C.mod_max` cycles are done before the physical
-        evolution time `_C.time_max` is reached the simulation stops
-        ahead of schedule.
+        should run. If `_C.mod_max` timesteps are done before the
+        physical evolution time `_C.time_max` is reached the simulation
+        stops ahead of schedule.
 
 -   `03` (`_C.freq_log`, `_C.freq_nir`, `_C.freq_ana`,
     `_C.freq_walltime`)
 
-    -   `_C.freq_log`: interval in units of timestep cycles at which the
-        NIRVANA log file `nirvana.log` (cf.
-        [Output data](3.3-Output-data#log-file)) and monitoring file
-        `nirvana.mon` (cf. [Output data](3.3-Output-data#monitoring-file))
-        is updated.
+    -   `_C.freq_log`: interval in units of timesteps at which the
+        NIRVANA log file `nirvana.log` and monitoring file `nirvana.mon`
+        (cf. [Output data](#Output-data)) is updated.
 
-    -   `_C.freq_nir`: interval in units of timestep cycles at which
-        NIRVANA snapshots (cf.
-        [Output data](3.3-Output-data#snapshot-files)) are produced.
+    -   `_C.freq_nir`: interval in units of timesteps at which NIRVANA
+        snapshots (cf. [Snapshot files](#snapshot-files)) are produced.
 
-    -   `_C.freq_ana`: interval in units of timestep cycles at which the
-        user interface function `analysisUser()` is called for
-        user-specific analysis tasks (cf.
-        [Data analysis](3.5-Data-analysis#runtime-data-analysis)).
+    -   `_C.freq_ana`: interval in units of timesteps at which the user
+        interface function `analysisUser()` is called for user-specific
+        analysis tasks (cf. [Data analysis](#data-analysis)).
 
     -   `_C.freq_walltime`: interval in seconds at which the special
-        snapshot `NIRLAST.#` (cf.
-        [Output data](3.3-Output-data#snapshot-files)) is renewed
-        overwriting existing older files.
+        snapshot `NIRLAST.#` (cf. [Snapshot files](#snapshot-files)) is
+        renewed by overwriting existing older files.
 
 **GEOMETRY**:
 
@@ -202,9 +185,8 @@ for the type of coordinate system.
 
     -   `_C.omega[0-2]`: components of the angular velocity vector of a
         rotating frame of reference with respect to the inertial frame
-        (cf. physics guide
-        [Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/A-PhysicsGuide),
-        [PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/PhysicsGuide.pdf)).
+        (cf. [physics
+        guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
         In case of cylindrical- or spherical geometry only the first
         component `_C.omega[0]` is accepted as relevant representing
         rotation around the geometric axis.
@@ -219,35 +201,34 @@ for the type of coordinate system.
 
 -   `01` (`_C.lo[0]`, `_C.up[0]`, `_C.dim[0]`)
 
-    -   `_C.lo[0]`,\_C.up\[0\]: lower,upper *x*-coordinate of the
+    -   `_C.lo[0]`,\_C.up\[0\]: lower,upper $x$-coordinate of the
         computational domain.
 
     -   `_C.dim[0]`: number of *base-level* grid points in
-        *x*-direction. `_C.dim[0]` must be an integral factor of 4 and
-        excludes ghost cells (which are automatically added).
+        $x$-direction. `_C.dim[0]` must be an integral factor of 4, and
+        excludes ghost cells which are automatically added by the code.
 
 -   `02` (`_C.lo[1]`, `_C.up[1]`, `_C.dim[1]`)
 
-    -   `_C.lo[1]`,`_C.up[1]`: lower,upper *y*-coordinate of the
-        computational domain. In case of spherical geometry (*y* ≡ *θ*)
-        `_C.lo[1]`,`_C.up[1]` have to be specified in units of *π*.
+    -   `_C.lo[1]`,`_C.up[1]`: lower,upper $y$-coordinate of the
+        computational domain. In case of spherical geometry
+        ($y\equiv \theta$) `_C.lo[1]`,`_C.up[1]` have to be specified in
+        units of $\pi$.
 
     -   `_C.dim[1]`: number of *base-level* grid points in
-        *y*-direction. `_C.dim[1]` must be a multiple factor of 4 and
-        excludes ghost cells (which are automatically added).
+        $y$-direction. `_C.dim[1]` must be a multiple factor of 4.
 
 -   `03` (`_C.lo[2]`, `_C.up[2]`, `_C.dim[2]`)
 
-    -   `_C.lo[2]`,`_C.up[2]`: lower,upper *z*-coordinate of the
+    -   `_C.lo[2]`,`_C.up[2]`: lower,upper $z$-coordinate of the
         computational domain. In case of cylindrical- or spherical
-        geometry (*z* ≡ *ϕ*) `_C.lo[2]`,`_C.up[2]` have to be specified
-        in units of *π*.
+        geometry ($z\equiv \phi$) `_C.lo[2]`,`_C.up[2]` have to be
+        specified in units of $\pi$.
 
     -   `_C.dim[2]`: number of *base-level* grid points in
-        *z*-direction. `_C.dim[2]` must be a multiple factor of 4 and
-        excludes ghost cells (which are automatically added). If
-        `_C.dim[2]`=0 the simulation is assumed 2D (e.g. axisymmetric in
-        case of cylindrical- or spherical coordinates).
+        $z$-direction. `_C.dim[2]` must be a multiple factor of 4. If
+        `_C.dim[2]`=0 the simulation is assumed 2D, i.e., axisymmetric
+        in case of cylindrical- or spherical coordinates.
 
 -   `04` (`_C.partitioning_type`, `_C.bnx`, `_C.bny`, `_C.bnz`)
 
@@ -255,16 +236,15 @@ for the type of coordinate system.
         MPI simulations. The value SFC stands for a Space-Filling-Curve
         domain decomposition and the value BLOCK for a rectangular
         domain decomposition. Both decomposition types work for unigrid
-        simulations. For AMR simulations a specified value BLOCK is
-        ignored and a SFC-decomposition is automatically applied.
+        simulations. For AMR simulations BLOCK is ignored and a
+        SFC-decomposition is automatically used instead.
 
     -   `_C.bnx`,`_C.bny`,`_C.bnz`: number of domain subdivisions in
-        *x*, *y*, *z*-direction in case \_C.partitioning\_type=BLOCK.
-        Numbers must be chosen such that the grid dimension of
-        subdomains is a multiple factor of 4 in each coordinate
-        direction. Moreover, the total number of subdomains must equal
-        the number of MPI threads, i.e.
-        `_C.bnx`×`_C.bny`×`_C.bnz`=number of MPI threads.
+        $x,y,z$-direction in case \_C.partitioning_type=BLOCK. Numbers
+        must be chosen such that the grid dimension of subdomains is a
+        multiple factor of 4 in each coordinate direction. Moreover, the
+        total number of subdomains must equal the number of MPI threads,
+        `_C.bnx`$\times$`_C.bny`$\times$`_C.bnz`=#threads.
 
 **BOUNDARY CONDITIONS**:
 
@@ -275,22 +255,21 @@ for the type of coordinate system.
 
     -   `_C.bc[0-5]` ({I,O,A,M,R,C,P,D,F,U}): selection of boundary
         conditions. The type of boundary condition at a physical domain
-        boundary is characterized by a single uppercase. Uppercases for
-        the domain boundaries are grouped in a 6-letters word with the
-        individual letters representing, from left to right, the
-        lower-*x* (`_C.bc[0]`), upper-*x* (`_C.bc[1]`), lower-*y*
-        (`_C.bc[2]`), upper-*y* (`_C.bc[3]`), lower-*z* (`_C.bc[4]`) and
-        upper-*z* (`_C.bc[5]`) domain boundary. Possible boundary
-        condition types are (*u*<sub>*i*/*o*</sub>: inner/outer-domain
-        boundary values of a variable *u*):
-
-        -   `I`: inflow – fluid is allowed to flow inwards through the
+        boundary is characterized by a single capital letter. It are
+        grouped in a 6-letters word with the individual letter
+        representing, from left to right, the lower-$x$ (`_C.bc[0]`),
+        upper-$x$ (`_C.bc[1]`), lower-$y$ (`_C.bc[2]`), upper-$y$
+        (`_C.bc[3]`), lower-$z$ (`_C.bc[4]`) and upper-$z$ (`_C.bc[5]`)
+        domain boundary. Possible boundary condition types are
+        ($u_{i/o}$: inner/outer-domain boundary values of a variable
+        $u$):
+
+        -   `I`: inflow -- fluid is allowed to flow inwards through the
             corresponding domain boundary but not to flow outwards. This
-            means for the perpendicular velocity component
-            *u*<sub>*o*</sub>=± sgn*u*<sub>*i*</sub>⋅*u*<sub>*i*</sub>
-            where the ’+’(’-’)-sign stands for a lower (upper) domain
-            boundary. Zero derivative,
-            *u*<sub>*o*</sub>=*u*<sub>*i*</sub>, is assumed for other
+            means For the perpendicular velocity component this means
+            $u_o=\pm\operatorname{sgn}u_i\cdot u_i$ where the
+            '+'('-')-sign is valid for a lower (upper) domain boundary.
+            Zero derivative, $u_o=u_i$, is assumed for other
             cell-centroid quantities except the total energy density.
             The total energy density is set consistently under the
             assumption of a zero derivative for the thermal energy
@@ -299,61 +278,56 @@ for the type of coordinate system.
             magnetic field component is obtained by a divergence-free
             extrapolation.
 
-        -   `O`: outflow – fluid is allowed to flow outwards through the
-            corresponding domain boundary but not to flow inwards. This
-            means for the perpendicular velocity component
-            *u*<sub>*o*</sub>=∓sgn *u*<sub>*i*</sub>⋅*u*<sub>*i*</sub>
-            where the ’+’(’-’)-sign stands for a lower (upper) domain
-            boundary. Other variables are likewise to I.
+        -   `O`: outflow -- fluid is allowed to flow outwards through
+            the corresponding domain boundary but not to flow inwards.
+            For the perpendicular velocity component this means
+            $u_o=\mp\operatorname{sgn}u_i\cdot u_i$ where the
+            '+'('-')-sign is valid for a lower (upper) domain boundary.
+            Other variables are set likewise to inflow I.
 
-        -   `M`: mirror symmetry – reflecting conditions,
-            *u*<sub>*o*</sub>=−*u*<sub>*i*</sub>, for the
-            perpendicular components of velocity and magnetic field and
-            zero derivative for the other variables.
+        -   `M`: mirror symmetry -- reflecting conditions, $u_o=-u_i$,
+            for the perpendicular components of velocity and magnetic
+            field. Zero derivative for the other variables.
 
-        -   `A`: antimirror symmetry – same as M for non-magnetic
+        -   `A`: anti-mirror symmetry -- same as M for non-magnetic
             variables. The magnetic field is forced to have dipole
             parity with respect to the domain boundary. This means
-            *u*<sub>*o*</sub>=−*u*<sub>*i*</sub> for the parallel
-            magnetic field components and
-            *u*<sub>*o*</sub>=*u*<sub>*i*</sub> for the perpendicular
-            magnetic field component.
+            $u_o=-u_i$ for the parallel magnetic field components and
+            $u_o=u_i$ for the perpendicular magnetic field component.
 
         -   `R`: reflection-on-axis (only of relevance in cylindrical
-            geometry for `_C.bc[2]` and in spherical geometry for
-            `_C.bc[2]`,`_C.bc[3]`) – reflecting conditions at the
-            geometric axis. Same as M except for the azimutal magnetic
-            field component for which
-            *u*<sub>*o*</sub>=−*u*<sub>*i*</sub>.
+            geometry at the $y$-lower boundary (at $R=0$) or in
+            spherical geometry at the $y$-lower/upper boundaries (at
+            $\theta=0,\pi$)) -- reflecting conditions at the geometric
+            axis. Same as M except for the azimutal magnetic field
+            component for which $u_o=-u_i$.
 
         -   `C`: reflection-at-center (only of relevance in spherical
-            geometry for `_C.bc[0]`) – reflecting conditions at the
-            coordinate origin. Same as M except for the non-radial
-            magnetic field components for which
-            *u*<sub>*o*</sub>=−*u*<sub>*i*</sub>.
-
-        -   `P`: periodicity – periodic conditions for all variables.
-
-        -   `D`: default – zero derivative in the non-magnetic variables
-            and perpendicular-to-boundary condition (pseudo-vacuum
-            condition) for the magnetic field.
-
-        -   `F`: free boundary (for cylindrical/spherical geometry
-            including the geometric axis) – ’natural’ boundary
-            conditions at the geometric axis. Boundary values are not
-            set explictely but are implicitly given by *π*-shifted
-            values at the geometric axis. Relevant in cylindrical
-            geometry for `_C.bc[2]` and in spherical geometry for
-            `_C.bc[0]`,`_C.bc[2]`,`_C.bc[3]`.
-
-            *Note 1: If the geometric axis is actually not part of the
-            domain setting F-type BC mimicks a pseudo-axis.*
-
-            *Note 2: There are restrictions in using F-type BC
-            concerning MPI and AMR.*
-
-        -   `U`: user-defined boundary conditions (cf.
-            [Defining BC](3.2-User-interfaces#defining-boundary-conditions)).
+            geometry at the $x$-lower boundary (at $r=0$)) -- reflecting
+            conditions at the coordinate center. Same as M except for
+            the non-radial magnetic field components for which
+            $u_o=-u_i$.
+
+        -   `P`: periodicity -- periodic conditions for all variables.
+
+        -   `D`: default -- zero derivative in the non-magnetic
+            variables and perpendicular-to-boundary condition
+            (pseudo-vacuum condition) for the magnetic field.
+
+        -   `F`: free boundary (only of relevance in cylindrical
+            geometry at the $y$-lower boundary (at $R=0$) or in
+            spherical geometry at the $y$-lower/upper boundaries (at
+            $\theta=0,\pi$) and $x$-lower boundary (at $r=0$)) --
+            'natural' boundary condition at the geometric axis. Boundary
+            values are not set explictely but are implicitly given by
+            $\pi$-shifted values. Note that when a `F`-type boundary
+            condition is used at a domain boundary not matching the
+            geometrical axis exactly a pseudo-axis is pretented. Note
+            also that restrictions apply in using F-type boundary
+            conditions in combination with MPI and AMR.
+
+        -   `U`: user-defined boundary conditions (cf. [Defining
+            boundary conditions](#defining-boundary-conditions)).
 
 **MESH REFINEMENT**:
 
@@ -366,26 +340,27 @@ for the type of coordinate system.
 
 -   `01` (`_C.imr`, `_C.amr`)
 
-    -   `_C.imr` (≤`MAXLEVEL`): maximal requested refinement level in a
-        user-defined initially refined mesh (cf.
-        [User-defined mesh refinement](3.2-User-interfaces#user-defined-mesh-refinement)).
+    -   `_C.imr` ($\le$`MAXLEVEL`): maximum refinement level for an
+        initially refined mesh by the user (cf. [User-defined initial
+        mesh refinement and refinement
+        control](#user-defined-initial-mesh-refinement-and-refinement-control)).
         `_C.imr` cannot be larger than the macro `MAXLEVEL` defined in
         the header file `nirvanaUser.h`.
 
-    -   `_C.amr` (≤`MAXLEVEL`): maximum allowed mesh refinement level in
-        a AMR simulations.
+    -   `_C.amr` ($\le$`MAXLEVEL`): allowed maximum mesh refinement
+        level in a AMR simulations. `_C.amr` cannot be larger than the
+        macro `MAXLEVEL` defined in the header file `nirvanaUser.h`.
 
 -   `02` (`_C.amr_eps[0-4]`)
 
-    -   `_C.amr_eps[0-4]` (typical range: \[0.1, 0.5\]): thresholds for
+    -   `_C.amr_eps[0-4]` (typical range: $[0.1,0.5]$): thresholds for
         the mass density (`_C.amr_eps[0]`), momentum density
         (`_C.amr_eps[1]`), energy density (`_C.amr_eps[2]`), magnetic
         field strength (`_C.amr_eps[3]`) and tracer variables
         (`_C.amr_eps[4]`) in the derivatives-based mesh refinement
-        criterion (cf.
-        [AMR](3.1-Code-basics#adaptive-mesh-refinement)). A zero or
-        negative value means that the respective component is disabled
-        for a mesh refinement check.
+        criterion (cf. [AMR](#adaptive-mesh-refinement)). A zero or
+        negative value indicates that the respective component is not
+        used in the mesh refinement check.
 
 -   `03` (`_C.rhoref`, `_C.vref`, `_C.eref`, `_C.bref`, `_C.Cref`)
 
@@ -415,42 +390,45 @@ for the type of coordinate system.
 
     -   `_C.amr_d1` (\[0,1\]): tuning parameter in the derivatives-based
         mesh refinement criterion controlling the transition between a
-        purely gradient-based (`_C.amr_d1=1`) and a purely
-        2nd-derivatives-based (`_C.amr_d1=0`) check.
+        purely gradient-based criterion (`_C.amr_d1=1`) and a purely
+        2nd-derivatives-based criterion (`_C.amr_d1=0`).
 
     -   `_C.amr_exp` (\[0,2\]): tuning parameter in the
         derivatives-based mesh refinement criterion controlling the
         power-law dependence of mesh refinement on the grid spacing. A
-        value `_C.amr_exp>1` (`_C.amr_exp<1`) would mean that mesh
-        refinement becomes progressively more difficult (easier) with
-        increasing refinement level compared to the standard linear
-        dependence (`_C.amr_exp=1`).
+        value `_C.amr_exp>1` (`_C.amr_exp<1`) means that mesh refinement
+        becomes progressively more difficult (easier) with increasing
+        refinement level compared to the standard linear dependence
+        (`_C.amr_exp=1`) (cf. [AMR](#adaptive-mesh-refinement)).
 
 -   `05` (`_C.amr_Jeans`, `_C.amr_exp`, `_C.amr_Field`)
 
-    -   `_C.amr_Jeans` (typical value: 0.2): threshold in the
+    -   `_C.amr_Jeans` (typical value: $0.2$): threshold in the
         Jeans-length-based mesh refinement criterion (cf.
-        [AMR](3.1-Code-basics#adaptive-mesh-refinement)). The value
-        `_C.amr_Jeans` defines the fraction of local Jeans length to be
-        resolved by at least one grid cell. A zero or negative value
-        means that the Jeans-length-based criterion is disabled.
+        [AMR](#adaptive-mesh-refinement)). The value `_C.amr_Jeans`
+        defines the fraction of local Jeans length to be resolved by one
+        grid cell. A zero or negative value means that the
+        Jeans-length-based criterion is disabled.
 
-    -   `_C.amr_dJeans` (≥0): tuning parameter for the
+    -   `_C.amr_dJeans` ($\ge 0$): tuning parameter for the
         Jeans-length-based mesh refinement criterion allowing a
         systematic reduction of the Jeans threshold with increasing
-        refinement level *l* according to the expression `_C.amr_Jeans`
-        -*l*\*`_C.amr_dJeans`, i.e., the local Jeans length becomes
-        gradually higher resolved with incresing *l*.
+        refinement level $l$ according to the expression `_C.amr_Jeans`
+        -$l$\*`_C.amr_dJeans`, i.e., the local Jeans length becomes
+        gradually higher resolved with incresing $l$. `_C.amr_dJeans`
+        must be positiv.
 
-        *Note: Tuning of `_C.amr_dJeans` must be such that the actual
-        Jeans threshold does not become too small of even negative.*
+        **Important**: `_C.amr_dJeans` must be chosen with care such
+        that the actual Jeans threshold never becomes too small or
+        negative at some refinement level. Otherwise it will trigger
+        mesh refinement up to the maximum refinement level allowed.
 
-    -   `_C.amr_Field` (typical value: 0.2): threshold in the
+    -   `_C.amr_Field` (typical value: $0.2$): threshold in the
         Field-length-based mesh refinement criterion (cf.
-        [AMR](3.1-Code-basics#adaptive-mesh-refinement)). The value
-        `_C.amr_Field` defines the fraction of the local Field length to
-        be resolved by at least one grid cell. A zero or negative value
-        means that the Field-length-based criterion is disabled.
+        [AMR](#adaptive-mesh-refinement)). The value `_C.amr_Field`
+        defines the fraction of the local Field length to be resolved by
+        one grid cell. A zero or negative value means that the
+        Field-length-based criterion is disabled.
 
 **USER-SPECIFIC PARAMETERS**:
 
@@ -462,10 +440,10 @@ for the type of coordinate system.
     05   0.00000e+00  0.00000e+00  0.00000e+00   >param[12-14] /
     06    0  0  0  0  0  0  0  0                 >flag[0-7] /
 
--   `01` (`_C.param[0-2]`)  
-    `02` (`_C.param[3-5]`)  
-    `03` (`_C.param[6-8]`)  
-    `04` (`_C.param[9-11]`)  
+-   `01` (`_C.param[0-2]`)\
+    `02` (`_C.param[3-5]`)\
+    `03` (`_C.param[6-8]`)\
+    `04` (`_C.param[9-11]`)\
     `05` (`_C.param[12-14]`)
 
     -   `_C.param[0-14]`: freely usable parameter of `double` type.
@@ -510,15 +488,15 @@ for the type of coordinate system.
         constrained-transport MHD solver.
 
         -   CT: flow-upwinded face-to-edge interpolation of the
-            face-centroid electric field fluxes obtained in the FV flux
+            face-centered electric field fluxes obtained in the FV flux
             computation part.
 
         -   CCT: 2D upwinding procedure based on a Hamilton-Jacobi-type
             numerical technique (equiv. to genuinely 2D-HLL in the
             current implementation).
 
-    -   `_C.mhd_courant` (typical value:  ≤ 0.5): CFL number in the
-        expression for the MHD timestep.
+    -   `_C.mhd_courant` (typical value: $\le 0.5$): CFL number in the
+        MHD timestep.
 
 -   `02` (`_C.viscosity_solver`, `_C.viscosity_courant`)
 
@@ -526,15 +504,14 @@ for the type of coordinate system.
         the fluid viscosity.
 
         -   STD: Runge-Kutta method as selected for the MHD integrator.
-            Coupling to the MHD integrator is unsplit.
+            Coupling to MHD integrator is unsplit.
 
         -   RKL: second-order, stabilized Runge-Kutta-Legendre method.
-            Coupling to the MHD integrator is via Strang-type splitting.
+            Coupling to MHD integrator is via Strang-type splitting.
 
-    -   `_C.viscosity_courant`: CFL-analog number in the expression for
-        the fluid viscosity timestep. A typical value in case of STD is
-         &lt;0.5. Values much larger than 1 are possible in the case of
-        RKL.
+    -   `_C.viscosity_courant`: CFL-like number in the fluid viscosity
+        timestep. A typical value in case of STD is $<0.5$. Values much
+        larger than 1 are possible in the case of RKL.
 
 -   `03` (`_C.diffusion_solver`, `_C.diffusion_courant`)
 
@@ -550,24 +527,26 @@ for the type of coordinate system.
 
 -   `06` (`_C.heatloss_max_change`)
 
-    -   `_C.heatloss_max_change` (typical value: ≤0.1): maximal
-        allowed relative change in temperature due to the heatloss
+    -   `_C.heatloss_max_change` (typical value: $\le 0.1$): allowed
+        maximal relative change in the temperature due to the heatloss
         source term.
 
 -   `07` (`_C.reactions_max_changeX`, `_C.reactions_max_changeT`)
 
-    -   `_C.reactions_max_changeX` (typical value: ≤0.1): maximal
-        relative change of species number densities (or total number
-        density) in the time integration of the chemo-thermal rate
-        equations.
+    -   `_C.reactions_max_changeX` (typical value: $\le 0.1$): allowed
+        maximal relative change of species number densities (or total
+        number density) in the time integration of the chemo-thermal
+        rate equations.
 
-        *Note: The macro `SXN` in module `solveNCCM.c` provides a linear
-        weight in applying the parameter `_C.reactions_max_change` to
-        individual number densities and to the total number density.*
+        **Important:** The macro `SXN` in module `solveNCCM.c` is a
+        linear weight applying `_C.reactions_max_change` to a combined
+        criterion based on individual number densities and the total
+        number density (`SXN`=1 uses individual number densities;
+        `SXN`=0 uses the total number density).
 
-    -   `_C.reactions_max_changeT` (typical value: ≤0.1): maximal
-        relative change of temperature in the time integration of the
-        chemo-thermal rate equations.
+    -   `_C.reactions_max_changeT` (typical value: $\le 0.1$): allowed
+        maximal relative change in the temperature in the time
+        integration of the chemo-thermal rate equations.
 
 -   `08` (`_C.dt0_reduce`)
 
@@ -592,137 +571,131 @@ for the type of coordinate system.
 
 -   `01` (`_C.mf`, `_C.permeability_rel`)
 
-    -   `_C.mf` ({Y,N}): option for activating magnetic fields. The
-        value Y enables the evolution of magnetic fields. The value N
-        refers to HD where the Lorentz force is dropped from the
-        momentum equation and the induction equation is not solved.
+    -   `_C.mf` ({Y,N}): option to include magnetic fields (MHD). Y
+        enables MHD. N refers to HD where the Lorentz force is dropped
+        from the momentum equation and the induction equation is not
+        solved.
 
     -   `_C.permeability_rel`: relative magnetic permeability
-        *μ*<sub>*r**e**l*</sub>
-        ( = *μ*/*μ*<sub>0</sub>, *μ*<sub>0</sub> = 4*π*⋅10<sup>−7</sup>V⋅m<sup>−1</sup>⋅A<sup>−1</sup>⋅s<sup>−1</sup>).
+        $\mu_{rel}$ ($=\mu/\mu_0,\,\mu_0=4\pi\cdot10^{-7} 
+        \mathtt{V}\cdot\mathtt{m}^{-1}\cdot\mathtt{A}^{-1}\cdot\mathtt{s}^{-1}$).
 
-        *Note: The Gaussian unit system can be mimicked by choosing a
-        value
-        `_C.permeability_rel`=10<sup>7</sup>
-        so that the magnetic permeability is *μ*=4*π*.*
+        **Important:** The Gaussian unit system can be mimicked by
+        choosing a value $\mathtt{\_C.permeability\_rel}=10^7$ so that
+        the magnetic permeability is $\mu=4\pi$.
 
 -   `02` (`_C.viscosity`, `_C.viscosity_coeff`)
 
-    -   `_C.viscosity` ({Y,N,U}): option for activating fluid viscosity.
-        The value Y enables the fluid viscosity term using a constant
-        *dynamic* viscosity coefficient as given by the parameter
-        `_C.viscosity_coeff` described next. The value U enables fluid
-        viscosity with a user-defined dynamic viscosity coefficient as
-        coded in the user interface `viscosityCoeffUser.c`. The value N
-        disables fluid viscosity.
+    -   `_C.viscosity` ({Y,N,U}): option for fluid viscosity. Y enables
+        the fluid viscosity term using a constant *dynamic* viscosity
+        coefficient as given by the parameter `_C.viscosity_coeff`
+        described below. U enables fluid viscosity with a user-defined
+        dynamic viscosity coefficient as coded in the user interface
+        function `viscosityCoeffUser()` (cf. [Fluid
+        viscosity](#fluid-viscosity)). N disables fluid viscosity.
 
-    -   `_C.viscosity_coeff`: standard constant *dynamic* viscosity
-        coefficient.
+    -   `_C.viscosity_coeff`: constant *dynamic* viscosity coefficient.
 
 -   `03` (`_C.diffusion`, `_C.diffusion_coeff`)
 
-    -   `_C.diffusion` ({Y,N,U}): option for activating magnetic
-        diffusion. The value Y enables the magnetic diffusion term using
-        a constant magnetic diffusion coefficient given by the parameter
-        `_C.diffusion_coeff` described next. The value U enables
-        magnetic diffusion with a user-defined magnetic diffusion
-        coefficient as coded in the user interface
-        `diffusionCoeffUser.c`. The value N disables magnetic
-        diffusion.
+    -   `_C.diffusion` ({Y,N,U}): option for Ohmic diffusion. Y enables
+        Ohmic diffusion using a constant magnetic diffusion coefficient
+        given by the parameter `_C.diffusion_coeff` described below. U
+        enables Ohmic diffusion with a user-defined magnetic diffusion
+        coefficient as coded in the user interface function
+        **diffusionCoeffUser()**. (cf. [Ohmic
+        diffusion](#ohmic-diffusion)). N disables Ohmic diffusion.
 
-    -   `_C.diffusion_coeff`: standard constant magnetic diffusion
-        coefficient.
+    -   `_C.diffusion_coeff`: constant magnetic diffusion coefficient.
 
 -   `04` (`_C.conduction`, `_C.conduction_coeff`,
     `_C.conduction_coeff_perp`, `_C.conduction_coeff_sat`)
 
-    -   `_C.conduction` ({Y,N,U,S}): option for activating thermal
-        conduction. In the presence of magnetic fields the thermal
-        conduction is anisotropic with distinct conduction coefficients
-        parallel and perpendicular to the magnetic field. In HD thermal
-        conduction is assumed isotropic. The value Y enables the thermal
-        conduction term using a constant conduction coefficient given by
-        the parameter `_C.conduction_coeff` (and
-        `_C.conduction_coeff_perp` in the anisotropic case) described
-        next. The value U enables thermal conduction with a user-defined
-        conduction coefficient as coded in the user interface
-        `conductionCoeffUser.c`. The value S refers to the standard
-        Spitzer conductivity model (cf. physics guide
-        [Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/A-PhysicsGuide),
-        [PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/PhysicsGuide.pdf)).
-        The value N disables thermal conduction.
-
-    -   `_C.conduction_coeff`: standard constant thermal conduction
-        coefficient (in the anisotropic case this is the parallel
-        coefficient).
-
-    -   `_C.conduction_coeff_perp`: constant thermal conduction
-        coefficient perpendicular to the magnetic field (meaningless in
-        isotropic conduction).
-
-    -   `_C.conduction_coeff_sat` (typical value: 0.3): parameter *Ψ* in
-        the saturation heat flux model of \[[CM77](3.2-User-interfaces#references)\] (cf. physics guide
-        [Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/C-NumericsGuide),
-        [PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/NumericsGuide.pdf)).
+    -   `_C.conduction` ({Y,N,U,S}): option for thermal conduction. In
+        the presence of magnetic fields the thermal conduction is
+        anisotropic with distinct conduction coefficients parallel and
+        perpendicular to the magnetic field. In HD thermal conduction is
+        assumed isotropic. Y enables thermal conduction using a constant
+        conduction coefficient given by the parameter
+        `_C.conduction_coeff` (and `_C.conduction_coeff_perp` in the
+        anisotropic case) described below. U enables thermal conduction
+        with a user-defined conduction coefficient as coded in the user
+        interface function **conductionCoeffUser()** (cf. [Thermal
+        conduction](#thermal-conduction)). S enables thermal conduction
+        using the standard Spitzer conductivity model predefined in
+        NIRVANA (cf. [physics
+        guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
+        N disables thermal conduction.
+
+    -   `_C.conduction_coeff`: constant thermal conduction coefficient
+        (in the anisotropic case this is the parallel coefficient).
+
+    -   `_C.conduction_coeff_perp`: thermal conduction coefficient
+        perpendicular to the magnetic field (meaningless in isotropic
+        conduction).
+
+    -   `_C.conduction_coeff_sat` (typical value: $0.3$): $\Psi$ is the
+        parameter in the saturation heat flux model of [@CM77] (cf.
+        [physics
+        guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
 
 -   `05` (`_C.APdiffusion`, `_C.APdiffusion_coeff`)
 
-    -   `_C.APdiffusion` ({Y,N,U}): option for activating ambipolar
-        diffusion. The value Y enables the ambipolar diffusion term
-        using a constant ambipolar diffusion coefficient given by the
-        parameter `_C.APdiffusion_coeff` described next. The value U
-        enables ambipolar diffusion with a user-defined anbipolar
-        diffusion coefficient as coded in the user interface
-        `APdiffusionCoeffUser.c`. The value N disables ambipolar
-        diffusion.
-
-    -   `_C.APdiffusion_coeff`: standard constant ambipolar diffusion
+    -   `_C.APdiffusion` ({Y,N,U}): option for ambipolar diffusion. Y
+        enables ambipolar diffusion using a constant ambipolar diffusion
+        coefficient given by the parameter `_C.APdiffusion_coeff`
+        described below. U enables ambipolar diffusion with a
+        user-defined anbipolar diffusion coefficient as coded in the
+        user interface function **APdiffusionCoeffUser()** (cf.
+        [Ambipolar
+        diffusion](#user-defined-coefficient-for-ambipolar-diffusion)).
+        N disables ambipolar diffusion.
+
+    -   `_C.APdiffusion_coeff`: constant ambipolar diffusion
         coefficient.
 
 -   `06` (`_C.energy`, `_C.energy_dual_sw`)
 
-    -   `_C.energy` ({Y,N,DUAL}): option for activation of the energy
-        equation. The value Y means that the total energy equation is
-        solved in time with the gasdynamics equations. For value N no
-        energy equation is solved. The value DUAL enables a dual energy
-        formalism where both the total energy equation and the thermal
-        energy equation is solved simultaneously. This allows a more
-        robust computation of the pressure (in view of possitivity) in
-        regions of high Mach number and/or low plasma-*β*. Actually, the
+    -   `_C.energy` ({Y,N,DUAL}): option for the energy equation. Y
+        means that a total energy equation is solved in step with
+        gasdynamics. For N no energy equation is solved. DUAL enables
+        the dual energy formalism where both the total energy equation
+        and the thermal energy equation is solved simultaneously. This
+        allows a more robust computation of the pressure in regions of
+        high Mach number and/or low plasma-$\beta$. In DUAL mode the
         thermal energy density (hence, the pressure) is computed from
         the solution of the thermal energy equation if the
         thermal-to-total energy density ratio falls below a certain
         threshold (next parameter). Otherwise the thermal energy density
-        is computed from the solution of the total energy equation as
-        usual.
-
-        *Note 1: Disabling the energy equation is not consistent with
-        the use of an adiabatic EOS.*
+        is computed from the solution of the total energy equation. Note
+        that the dual energy formalism violates the energy conservation
+        principle since syncronization between the total- and thermal
+        energy density is necessary.
 
-        *Note 2: The dual energy formalism violates the energy
-        conservation principle since syncronization steps between total
-        and thermal energy densities are necessary to maintain
-        consistency.*
+        **Important:** Disabling the energy equation is not compatible
+        with the use of an adiabatic EOS, for instance.
 
-    -   `_C.energy_dual_sw` (\[0, 1\]; typical value: 0.01): threshold
+    -   `_C.energy_dual_sw` ($[0,1]$; typical value: $0.01$): threshold
         value for the thermal-to-total energy density ratio in the dual
         energy formalism (when `_C.energy`=DUAL).
 
 -   `07` (`_C.force`)
 
-    -   `_C.force` ({Y,N}): option for activating a body force term. The
-        value Y enables a body force in the momentum equation. The body
-        force has to be defined by the user in the user interface
-        `forceUser.c`. The value N disables body forces.
+    -   `_C.force` ({Y,N}): option for an external body force term. Y
+        enables a body force in the momentum equation. The body force
+        has to be defined by the user in the user interface function
+        `forceUser()` (cf. [User-defined body
+        force](#user-defined-body-force)). N disables body forces.
 
 -   `08` (`_C.gravity`)
 
-    -   `_C.gravity` ({Y,N}): option for activating selfgravity. The
-        value Y enables selfgravity. A Poisson equation is solved for
-        the gravitational potential (cf. physics guide
-        [Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/C-NumericsGuide),
-        [PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/NumericsGuide.pdf)).
-        The value N disables selfgravity.
+    -   `_C.gravity` ({Y,N,S}): option for selfgravity. Y enables
+        selfgravity. S enables selfgravity using the screening mass
+        approach to solve for Dirichlet boundary conditions. S works
+        only for cylindrical/spherical geometry with $2\pi$ periodicity
+        in azimuthal direction (cf. [physics
+        guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
+        N disables selfgravity.
 
 -   `09` (`_C.eos`,`_C.gamma`,`_C.polytropic_constant`,`_C.temperature`)
 
@@ -730,33 +703,33 @@ for the type of coordinate system.
         state. Possible values are:
 
         -   ADIA: adiabatic EOS. This choice requires specificiation of
-            the adiabatic index (parameter: `_C.gamma`).
+            the adiabatic index `_C.gamma`.
 
         -   ISO: isothermal EOS. This choice requires specification of
-            the temperature (parameter: `_C.temperature`).
+            the temperature `_C.temperature`.
 
         -   POLY: polytropic EOS. This choice requires specification of
-            the polytropic exponent (parameter: `_C.gamma`) and
-            polytropic constant (parameter: `_C.polytropic_constant`).
+            the polytropic exponent `_C.gamma` and polytropic constant
+            `_C.polytropic_constant`.
 
         -   USER: user-defined analytic EOS. This choice requires the
-            specification of special macros (cf.
-            [User-defined EOS](3.2-User-interfaces#user-defined-equation-of-state))
+            specification of special macros (cf. [User-defined
+            EOS](#user-defined-equation-of-state))
 
         -   TAB: user-defined tabulated EOS. This choice requires the
-            generation of look-up tables (cf.
-            [User-defined EOS](3.2-User-interfaces#user-defined-equation-of-state)).
+            generation of look-up tables (cf. [User-defined
+            EOS](#user-defined-equation-of-state)).
 
-        *Note: The choice of an isothermal EOS or polytropic EOS is not
-        compatible with solving for an energy equation.*
+        **Important:** The choice of an isothermal EOS or polytropic EOS
+        is not compatible with solving for an energy equation.
 
     -   `_C.gamma`: adiabatic index (polytropic exponent) in the
-        adiabatic (polytropic) EOS.
+        adiabatic (polytropic) EOS. Note that in a multi-species gas the
+        adiabatic index is not a freely selectable parameter but is
+        approximated by the expression
 
-        *Note: In a multi-species gas the adiabatic index is not a
-        freely selectable parameter but is approximated by the
-        expression*  
-        *γ*=\[5(*n*<sub>*H*</sub>+*n*<sub>*He*</sub>+*n*<sub>*e*</sub>)+7*n*<sub>*H*<sub>2</sub></sub>\]/\[3(*n*<sub>*H*</sub>+*n*<sub>*He*</sub>+*n*<sub>*e*</sub>)+5*n*<sub>*H*<sub>2</sub></sub>\]
+        $\gamma=\left[5(n_\mathrm{H}+n_\mathrm{He}+n_\mathrm{e})+7n_{\mathrm{H}_2}\right]
+        /\left[3(n_\mathrm{H}+n_\mathrm{He}+n_\mathrm{e})+5n_{\mathrm{H}_2}\right]$.
 
     -   `_C.polytropic_constant`: polytropic constant in the polytropic
         EOS.
@@ -765,18 +738,17 @@ for the type of coordinate system.
 
 -   `10` (`_C.heatloss`)
 
-    -   `_C.heatloss` ({N,ISM,U}): option for activating an external
-        heatloss term (cooling/heating function) in the energy equation.
-        The value N disables the heatloss term. The value ISM enables
-        the heatloss term applying a predefined interstellar medium
-        cooling/heating function (cf. physics guide
-        [Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/A-PhysicsGuide),
-        [PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/PhysicsGuide.pdf)).
-        The value U enables the heatloss term assuming a user-defined
+    -   `_C.heatloss` ({N,ISM,U}): option for an external heatloss
+        source term (cooling/heating function) in the energy equation. Y
+        enables the heatloss term. N disables it. ISM enables the
+        heatloss term and applies a predefined interstellar medium
+        cooling/heating function (cf. [physics
+        guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
+        U enables the heatloss term assuming a user-defined
         cooling/heating function coded in the interfaces
-        `sourceCoolingUser.c` for cooling and `sourceHeatingUser.c`
-        for heating (cf.
-        [User-defined cooling/heating function](3.2-User-interfaces#user-defined-cooling-and-heating-function)).
+        **sourceCoolingUser()** for cooling and **sourceHeatingUser()**
+        for heating (cf. [User-defined
+        cooling/heating](#user-defined-cooling-and-heating-function)).
 
 -   `11` (`_C.tracer`)
 
@@ -785,108 +757,125 @@ for the type of coordinate system.
 
 -   `12` (`_C.testfields`)
 
-    -   `_C.testfields`: number of testfield equations associated with
-        the code functionality TESTFIELDS.  
+    -   `_C.testfields`: number of testfields associated with the code
+        functionality TESTFIELDS.
 
--   `13` (`NN`, `NN`, `_C.energy_reactions`, `_C.mean_molecular_weight`,
+-   `13` ( NN, NN, `_C.energy_reactions`, `_C.mean_molecular_weight`,
     `_C.mean_ionization`)
 
-    -   NN ({Y,N}): option for activating the multi-species framework.
-        The value Y enables the multi-species framework. In the
-        multi-species framework an advection equation is solved for each
-        specie defined in the `SPECIES` section of the user interface
-        `NCCM.par`. The value N refers to a standard single-component
-        gas.
-
-    -   NN ({Y,N}): option for activating a chemical reaction network
-        between species. The value Y enables chemical reactions if the
-        multi-species framework is enabled. An advection-reaction
-        equation is solved for each specie defined in the `SPECIES`
-        section of the user interface `NCCM.par` subject to chemical
-        reactions as specified in the `REACTIONS` section of the same
-        file.
-
-    -   `_C.energy_reactions` ({Y,N}): option for activating the thermal
-        processes implemented in the NCCM. The value Y enables thermal
-        processes if the multi-species framework and chemical reaction
-        network are both enabled. Thermal processes have to be selected
-        in the `THERMAL PROCESSES` section of `NCCM.par`.
+    -   NN ({Y,N}): option for a multi-species framework. Y enables a
+        multi-species framework. In the multi-species framework an
+        advection equation is solved for each gas specie. Species are
+        defined in the `SPECIES` section in the parameter file
+        `NCCM.par`. N refers to using a canonical single-species gas.
+
+    -   NN ({Y,N}): option for a chemical reaction network. Y enables
+        chemical reactions between species if a multi-species framework
+        is enabled. In this case an advection-reaction equation is
+        solved for each gas specie subject to chemical reactions.
+        Chemical reactions are specified in the `REACTIONS` section in
+        `NCCM.par`.
+
+    -   `_C.energy_reactions` ({Y,N}): option for chemo-thermal
+        processes implemented in the NCCM. Y enables thermal processes
+        if both a multi-species framework and chemical reaction network
+        is enabled. Thermal processes are selected in the
+        `THERMAL PROCESSES` section of `NCCM.par`. N disables
+        chemo-thermal processes.
 
     -   `_C.mean_molecular_weight`: mean molecular weight.
 
-        *Note: In a multi-species gas the mean molecular weight is not a
-        freely selectable parameter but computed selfconsistently from
-        the species composition.*
+        **Important:** In a multi-species framework the mean molecular
+        weight is not a freely selectable parameter but it is computed
+        selfconsistently from the gas composition.
 
     -   `_C.mean_ionization`: mean ionization fraction.
 
-        *Note: In a multi-species gas the ionization fraction is not a
-        freely selectable parameter but computed selfconsistently from
-        the ionisation structure.*
-
-## Defining initial conditions
-
-The user interface `configUser.c` serves to define IC for a problem.
-Defining IC means assigning values to primary physical variables,
-{𝜚, **m**, *e*, **B**, *C*<sub>*c*</sub>, *n*<sub>*s*</sub>}, on the
-mesh. However, only problem-relevant variables have to be defined.
-Primary variables which are irrelevant for the problem need not to be
-assigned. For example, in simulations with isothermal EOS the total
-energy density *e* is irrelevant since no energy equation is solved.
-Assignment of the magnetic field in HD simulations (when `_C.mf`=0) is
-redundant and even not allowed because arrays for the magnetic field
-components are not allocated in the HD case. Likewise, arrays for
-*C*<sub>*c*</sub> and *n*<sub>*s*</sub> are not declared if
-`_C.tracer`=0 respective `_C.species`=0. Furthermore, derived variables
-like the temperature must not be assigned here. If the TESTFIELDS
-infrastructure is used testfield fluctuation variables,
-**b**<sub>*t*</sub>, *t*=0,*N*<sub>*t*</sub>−1, are to be
-considered primary.
-
-*Note: In case magnetic field splitting is enabled
-(cf. [Functionality overview](3.1-Code-basics#functionality-overview))
-the variable* **B** *represents the residual magnetic field which
-must be initialized instead of the total field. Likewise, *e* represents
-the residual total energy density.*
-
-The mesh is represented by the master mesh pointer `gm` which is the
-only argument passed to function `configUser()`. The problem-relevant
-variables have to be assigned for each superblock `g` in `gm`. Recall
-in this context the explanations about the [mesh data
-structure](3.1-Code-basics\#mesh-data-structure).
-
-There are two types of primary variables: cell-averaged variables and
+        **Important:** In a multi-species framework the ionization
+        fraction is not a freely selectable parameter but is computed
+        selfconsistently from the ionisation degree of gas species.
+
+### Defining initial conditions
+
+The module `configUser.c` serves as user interface for the definition of
+initial conditions (IC). To do this the user must program the function
+`configUser()` in the module. Defining IC means assigning values to
+primary physical variables on the mesh, that are
+$\{\varrho,\mathbf{m},e,\mathbf{B},C_\mc,n_\ms\}$ -- gas density,
+momentum,total energy density, magnetic field, tracer
+($c=0,{\tt\_C.tracer}-1$), species number densities
+($s=0,{\tt\_C.species}-1$). The parameter `_C.tracer` (`_C.species`)
+denotes the number of tracer variables (species). However, those primary
+variables which are irrelevant for the problem under study need not to
+be assigned. For example, in simulations with isothermal EOS the total
+energy density $e$ is irrelevant since no energy equation is solved. In
+HD simulations, the assignment of the magnetic field is redundant, of
+course, and would lead to an error because arrays for the magnetic field
+components are not allocated in the HD case. Likewise, arrays for tracer
+and species, $C_\mc$ and $n_\ms$, are not declared unless the parameters
+`_C.tracer`$>0$ respective ${\tt\_C.species}>0$. If the TESTFIELDS
+infrastructure is used (`_C.testfields`$>0$) testfields fluctuation
+variables, $\mathbf{b}_\mt,\,\mt=0,{\tt\_C.testfields}-1$, are also
+considered as primary variables and must be assigned by the user. The
+parameter `_C.testfields` denotes the number of testfields.
+
+**Important:** In case magnetic field splitting is enabled (cf. [Code
+features](#code-features)) the variable $\mathbf{B}$ represents the
+residual magnetic field, i.e., total field minus background field. In
+addition, $e$ represents the residual total energy density, i.e, has a
+magnetic energy density contribution solely from the residual magnetic
+field.
+
+The function
+
+      configUser(gm)
+
+is passed the master mesh pointer `gm` which contains all the mesh and
+data structures. In serial simulations `gm` represents the global mesh.
+In MPI simulations, `gm` only represents the mesh partition on a thread.
+Primary variables have to be assigned for each mesh superblock `g` in
+`gm`. Looping over the mesh is looping over all refinement levels,
+il=0,`_C.level`, and all its superblocks:
+
+      for(il=0; il<=_C.level; il++)
+        for(g=gm[il]; (g); g=g->next)
+        {
+          ...
+        }
+
+Please recall in this context the introduction to NIRVANA's [Mesh data
+structure](#mesh-data-structure).
+
+There are two types of mesh variables: cell-averaged variables and
 face-averaged variables. Cell-averaged variables are
-𝜚, *m*<sub>*x*</sub>, *m*<sub>*y*</sub>, *m*<sub>*z*</sub>, *Φ*, *n*<sub>*s*</sub>, *C*<sub>*c*</sub>
-and it are represented by cell-centroid coordinates. Face-averaged
-variables are the magnetic field components
-*B*<sub>*x*</sub>, *B*<sub>*y*</sub>, *B*<sub>*z*</sub> and it are
-represented by face-centroid cell coordinates.
-
-A superblock `g` contains ghost cells added to its active region.
-However, in the definition of IC ghost cells need not to be assigned by
-the user. It is sufficient to assign active cells, i.e. those which are
-integrated. Active cells of a superblock start at lower grid index
-(`g->ixs`,`g->iys`,`g->izs`). The upper index of active cells depends on
-the type of variable and is given in the following table together with
-the representative location of the various variables in a cell:
-
-| variable                                                                                              | upper index                           | cell location                         |
+$\varrho,m_x,m_y,m_z,n_\ms,C_\mc$. It are described by cell-centroid
+coordinates. Face-averaged variables are the magnetic field components
+$B_x,B_y,B_z$ (and testfields). It are described by cell face-centroid
+coordinates.
+
+Although a superblock `g` contains ghost cells (taking care of boundary
+values) it is sufficient to only assign cells of the active region.
+Active cells of a superblock start at lower grid index
+(`g->ixs`,`g->iys`,`g->izs`). The upper index depends on the type of
+variable, whether cell-averaged or face-averaged, and is given in the
+following table. The table also compiles the coordinates of the
+variables in a cell:
+
+| variable                                                                                              | upper index                           | cell coordinate location              |
 |:------------------------------------------------------------------------------------------------------|:--------------------------------------|:--------------------------------------|
 | 𝜚, *m*<sub>*x*</sub>, *m*<sub>*y*</sub>, *m*<sub>*z*</sub>, *Φ*, *n*<sub>*s*</sub>, *C*<sub>*c*</sub> | (`g->ixe`,`g->iye`,`g->ize`)          | (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`) |
 | *B*<sub>*x*</sub>                                                                                     | (`g->ixe+1,g->iye,g->ize`)            | (`g->x[ix]`,`g->yc[iy]`,`g->zc[iz]`)  |
 | *B*<sub>*y*</sub>                                                                                     | (`g->ixe`,`g->iye`+1,`g->ize`)        | (`g->xf[ix]`,`g->y[iy]`,`g->zc[iz]`)  |
 | *B*<sub>*z*</sub>                                                                                     | (`g->ixe`,`g->iye`,`g->ize`+`_C.idz`) | (`g->xf[ix]`,`g->yc[iy]`,`g->z[iz]`)  |
 
-*Note 1: The staggered, face-averaged vector components of the magnetic
-field have an upper index increased by one in the respective coordinate
-direction.*
-
-*Note 2: The code parameter `_C.idz` indicates whether a problem is 2D
+Note that the face-averaged vector components of the magnetic field (and
+likewise the testfields components) have an upper index increased by one
+in the respective coordinate direction compared to the cell-averaged
+variables. The code parameter `_C.idz` indicates whether a problem is 2D
 (`_C.idz`=0; axissymmetry in cylindrical/spherical geometry) or 3D
-(`_C.idz`=1).*
+(`_C.idz`=1).
 
-Looping over the active grid cells of a superblock looks like this:
+Here is an example of looping over the active region of a superblock:
 
       /* ASSIGNMENT OF CELL-AVERAGED VARIABLES */
       for(iz=g->izs; iz<=g->ize; iz++)
@@ -895,6 +884,14 @@ Looping over the active grid cells of a superblock looks like this:
           { 
             g->rho[iz][iy][ix]=...
             ...
+
+            if(_C.tracer) /* SET TRACER VARIABLES */
+              for(ic=0; ic<_C.tracer; ic++)
+                g->C[ic][iz][iy][ix]=...
+
+            if(_C.species) /* SET SPECIES */
+              for(is=0; is<_C.species; is++)
+                g->nX[is][iz][iy][ix]=...
           }
 
       /* ASSIGNMENT OF FACE-AVERAGED MAGNETIC FIELD COMPONENTS */
@@ -916,48 +913,50 @@ Looping over the active grid cells of a superblock looks like this:
           for(ix=g->ixs; ix<=g->ixe; ix++)
             g->bz[iz][iy][ix]=...
 
-Since tracer variables, *C*<sub>*c*</sub>, are cell-averaged quantities
-like the mass density it can be assigned at the same place. The tracer
-array index, `ic`, runs from 0 to `_C.tracer`-1 where the code parameter
-`_C.tracer` gives the number of tracers *N*<sub>*c*</sub>. Tracer
-variables are dimensionless and are usually defined in the range
-\[0,1\].
+The tracer array index, `ic`, runs from 0 to `_C.tracer`-1 whereas the
+species array index, `is`, runs from 0 to `_C.species`-1. Tracer
+variables are assumed dimensionless and it are usually defined in the
+range $[0,1]$. The mesh of the superblock is indexed by
+(`ix`,`iy`,`iz`).
 
-Similary, species number densities, *n*<sub>*s*</sub>, are cell-averaged
-quantities. The species array index, `is`, runs from 0 to `_C.species`-1
-where code parameter `_C.species` gives the number of species
-*N*<sub>*s*</sub>. The species number densities must be consistently
-defined along with the mass density expressed by the constraint
+**Important:** The species number densities must be consistently defined
+along with the mass density because of the conservation condition
 
-∑<sub>*s*</sub>(*μ*<sub>*s*</sub> + *μ*<sub>*e*</sub>*q*<sub>*s*</sub>)*n*<sub>*s*</sub> = 𝜚/*m*<sub>*u*</sub>,  *s* = 0, *N*<sub>*s*</sub> − 1
+$$\sum_s(\mu_s+\mu_e q_s)n_s = \varrho/m_u,\quad s=0,{\tt \_C.species}-1$$
 
-where *μ*<sub>*s*</sub>, *q*<sub>*s*</sub> and *μ*<sub>*e*</sub> are the
-molecular weight and charge of the s-th specie and the electron weight,
-respectively. In principle, number densities can be assigned
-individually on this premise. Alternatively, and much simpler, number
-densities can be computed in compliance with the user-specified mass
-fractions in the parameter file `NCCM.par` by calling the function
+where $\mu_s$, $q_s$ and $\mu_\mathrm{e}$ are the molecular weight,
+charge of the s-th specie and the electron weight, respectively. In
+principle, number densities can be assigned individually in this way.
+However, it may be simpler to specify the species mass fractions in the
+parameter file `NCCM.par` and then to call the function
 
       CHinitSpecies(gm);
 
-in the *function body* of `configUser.c` after the mass density has been
-assigned. Consistency between the species number densities and mass
-density is thereby automatically guaranteed.
+in the *body* (outside of `g`-loops) of `configUser()` after the mass
+density has been assigned. It guarantees consistency between the species
+number densities and mass density.
 
-A frequently encounterd problem in the definition of IC for the magnetic
-field is to find a divergence-free formulation on the mesh consistent
-with the Constrained-Transport integral approach. This may be a
-non-trivial task for non-uniform fields, in particular on an initially
-refined mesh. There are two clean approaches to achieve a
+**On defining a divergence-free magnetic field.** A frequently
+encounterd problem in the definition of IC for the magnetic field is to
+fulfill the divergence contraint ($\nabla\cdot\mathbf{B}=0$) on the mesh
+in a Finite-Volume sense. This may be a non-trivial task, in particular,
+for non-uniform fields. There are two clean approaches to achieve a
 divergence-free setup:
 
-\(1\) computation of cell-face-averaged magnetic field components via
-exact integration of the analytical expressions. The so discretized
-field is per se cell-wise divergence-free. For a grid cell
-(`ix`,`iy`,`iz`) on superblock `g` the discritized components read
-
-![ic_b](uploads/411b3e817fec94afa0a2861772697e71/ic_b.png)  
-The numerical expressions for the cell face contents are:
+**(I)** Explicit computation of cell-face-averaged magnetic field
+components by exact integration of given analytical expressions. For a
+grid cell (,,) on superblock `g` the discretized components would then
+read 
+$$\mathtt{g->bx[iz][iy][ix]}=\frac{1}{\d \mathcal{A}_x}\int
+\limits_\mathtt{g->y[iy]}^\mathtt{g->y[iy+1]}
+\int\limits_\mathtt{g->z[iz]}^\mathtt{g->z[iz+1]} B_x(\mathtt{g->x[ix]},y,z)h_yh_zdydz$$
+$$\mathtt{g->by[iz][iy][ix]}=\frac{1}{\d \mathcal{A}_y}\int
+\limits_\mathtt{g->x[ix]}^\mathtt{g->x[ix+1]}
+\int\limits_\mathtt{g->z[iz]}^\mathtt{g->z[iz+1]} B_y(x,\mathtt{g->y[iy]},z)h_zdxdz$$
+$$\mathtt{g->bz[iz][iy][ix]}=\frac{1}{\d \mathcal{A}_z}\int
+\limits_\mathtt{g->x[ix]}^\mathtt{g->x[ix+1]}
+\int\limits_\mathtt{g->y[iy]}^\mathtt{g->y[iy+1]} B_z(x,y,\mathtt{g->z[iz]})h_ydxdy$$
+where the numerical expressions for the cell faces are:
 
 | cell face            | expression                               |
 |:---------------------|:-----------------------------------------|
@@ -965,68 +964,83 @@ The numerical expressions for the cell face contents are:
 | *δ*𝒜<sub>*y*</sub> | `g->dax[ix]*g->hzyh[iy]*g->dz`           |
 | *δ*𝒜<sub>*z*</sub> | `g->dax[ix]*g->dy`                       |
 
-\(2\) use of the magnetic vector potential **A**, if known, with
-**B** = ∇ × **A**. When discretized in integral form this gives for the
-face-averaged magnetic field components
+This way of approach is in most cases not feasible because the 2D
+integration over cell faces turns out to be too complicated.
+
+**(II)** Use of the magnetic vector potential, $\mathbf{A}$, if known.
+The face-averaged magnetic field components are then obtained from the
+integral form of $\mathbf{B}=\nabla\times \mathbf{A}$:
+
+`g->bx[iz][iy][ix]`$=\left[h_y\Delta_\miy (h_{zy}\hat{A}_z)
+-h_y\Delta_\miz \hat{A}_y\right]/\d \mathcal{A}_x$
+
+`g->by[iz][iy][ix]`$=\left[\Delta_\miz \hat{A}_x
+-h_{zy}\Delta_\mix (h_y\hat{A}_z)\right]/\d \mathcal{A}_y$
 
-`g->bx[iz][iy][ix]` = \[*h*<sub>*y*</sub>*Δ*<sub>`iy`</sub>(*h*<sub>*zy*</sub>*Â*<sub>*z*</sub>)−*h*<sub>*y*</sub>*Δ*<sub>`iz`</sub>*Â*<sub>*y*</sub>\]/*δ*𝒜<sub>*x*</sub>
+`g->bz[iz][iy][ix]`$=\left[\Delta_\mix (h_y\hat{A}_y)
+-\Delta_\miy \hat{A}_x\right]/\d \mathcal{A}_z$
 
-`g->by[iz][iy][ix]` = \[*Δ*<sub>`iz`</sub>*Â*<sub>*x*</sub>−*h*<sub>*zy*</sub>*Δ*<sub>`ix`</sub>(*h*<sub>*y*</sub>*Â*<sub>*z*</sub>)\]/*δ*𝒜<sub>*y*</sub>
+where the difference operators $\Delta_\mix$, $\Delta_\miy$ and
+$\Delta_\miz$ are given by
 
-`g->bz[iz][iy][ix]` = \[*Δ*<sub>`ix`</sub>(*h*<sub>*y*</sub>*Â*<sub>*y*</sub>)−*Δ*<sub>`iy`</sub>*Â*<sub>*x*</sub>\]/*δ*𝒜<sub>*z*</sub>
+$\Delta_{\mix}X=X(\mathtt{ix+1,iy,iz})-X(\mathtt{ix,iy,iz})$,
 
-where (*Â*<sub>*x*</sub>, *Â*<sub>*y*</sub>, *Â*<sub>*z*</sub>) denote
-the path integrals
-![ic_A](uploads/d7670a3203d0f3a1e2a9a759e69e2e77/ic_A.png)  
-over corresponding cell edges. The difference operator
-*Δ*<sub>`iy`</sub>, for instance, is given by
+$\Delta_{\miy}X=X(\mathtt{ix,iy+1,iz})-X(\mathtt{ix,iy,iz})$,
 
-*Δ*<sub>`iy`</sub>*Â*<sub>*x*</sub> = *Â*<sub>*x*</sub>(`ix`,`iy`+1,`iz`) − *Â*<sub>*x*</sub>(`ix`,`iy`,`iz`)
+$\Delta_{\miz}X=X(\mathtt{ix,iy,iz+1})-X(\mathtt{ix,iy,iz})$
 
-and similarly for others. If the *Â*’s are *unambiguously* computed (or
-approximated) on coinciding cell edges in the grid hierarchy, the
-magnetic field has a consistent, cell-wise divergence-free mesh
-representation.
+and $(\hat{A}_x,\hat{A}_y,\hat{A}_z)$ are the cell-edge integrals
 
-*Note 3: The code can not check for the consistency of user-defined
-magnetic field configurations. If inconsistencies exist initially,
-especially on refined meshes – even not obvious – it may lead to
-magnetic monopole generation during run time.*
+$$\hat{A}_x(\mathtt{ix,iy,iz})=\int\limits_\mathtt{g->x[ix]}^\mathtt{g->x[ix+1]}
+A_x(x,\mathtt{g->y[iy]},\mathtt{g->z[iz]})dx$$
 
-If the TESTFIELDS infrastructure is enabled and IC are needed for
-testfield fluctuation variables it must be handled in the same way as
-the magnetic field. Moreover, the same issues concerning a
-divergence-free mesh representation appear as for the magnetic field.
-The array index for the testfield fluctuation variables runs from 0 to
-`_C.testfields`-1 where the code parameter `_C.testfields` defines the
-number of testfields *N*<sub>*t*</sub>.
+$$\hat{A}_y(\mathtt{ix,iy,iz})=\int\limits_\mathtt{g->y[iy]}^\mathtt{g->y[iy+1]}
+A_y(\mathtt{g->x[ix]},y,\mathtt{g->z[iz]})dy$$
 
-The setup of IC in MPI parallel simulations is more complicated than in
-serial simulations because of mesh partitioning. The mesh master pointer
-`gm` passed to the IC user interface `configUser.c` only represents the
-portion of mesh distributed to the calling MPI thread. The user has to
-take into account here the spatial locations of superblocks when
-assigning physical data to the mesh.
+$$\hat{A}_z(\mathtt{ix,iy,iz})=\int\limits_\mathtt{g->z[iz]}^\mathtt{g->z[iz+1]}
+A_z(\mathtt{g->x[ix]},\mathtt{g->y[iy]},z)dz$$
 
-Many examples of how to define IC can be found in the testproblems
-collection located in subdirectories under `/nirvana/testproblems/`.
-Here, the IC for the Orszag-Tang problem (example 1) and a shock-cloud
-interaction problem (example 2) are presented as examples. Both are MHD
-problems.
+The $\hat{A}$'s are usually easier to compute than the face-averaged
+magnetic field components directly. Even if the $\hat{A}$'s are
+approximated but otherwise set *unambiguously* on cell edges in the grid
+hierarchy, the magnetic field derived from the above discrete curl of
+vector potential has a divergence-free mesh representation.
 
-**Example 1** (taken from `/nirvana/testproblems/MHD/problem2`; cf. \[[Zie04](3.2-User-interfaces#references)\])  
+Note that NIRVANA can not check for the consistency of a user-defined
+magnetic field configuration. If inconsistencies exist initially,
+especially on a refined mesh -- even not obvious -- it may lead to
+magnetic monopole generation during run time.
 
-IC for the Orszag-Tang problem simulated in a doubly-periodic square
-domain of length *L* (given by `_C.up[0]-_C.lo[0]`):
+When the TESTFIELDS infrastructure is enabled and IC are needed for
+testfields fluctuation variables it must be handled in the same way as
+the magnetic field. In particular, the same issue concerning a
+divergence-free representation on the mesh applies.
 
-𝜚 = 𝜚<sub>0</sub>  
-**m** = 𝜚<sub>0</sub>(−sin(2*π* *y*/*L*),sin(2*π* *x*/*L*),0)  
-**B** = *B*<sub>0</sub>(−sin(2*π* *y*/*L*),sin(4*π* *x*/*L*),0)  
-*e* = *p*<sub>0</sub>/(*γ*−1) + **m**<sup>2</sup>/(2𝜚) + **B**<sup>2</sup>/(2*μ*)
+**Important:** For MPI simulations in case the IC require global mesh
+operations parallel coding using MPI functionality is needed.
 
-where 𝜚<sub>0</sub>=2.77778, *p*<sub>0</sub>=5/3 and
-*B*<sub>0</sub>=1. The adiabatic index *γ*=5/3 and the magnetic
-permeability *μ*=1.
+Many examples of how to define IC can be found in NIRVANA's collection
+of testproblems. It are located in subdirectories below
+`/nirvana/testproblems/`. Here, the definition of IC are shown examplary
+for two problems: the Orszag-Tang vortex problem (example 1) and a
+shock-cloud interaction problem (example 2).
+
+**Example 1** (taken from `/nirvana/testproblems/MHD/problem2`; cf.
+[@Z04])
+
+IC for the Orszag-Tang vortex problem simulated in a doubly-periodic
+square domain of length $L$ (given by `_C.up[0]-_C.lo[0]`):
+
+$\varrho = \varrho_0$
+
+$\mathbf{m}=\varrho_0 \left(-\sin (2\pi\,\!y/L),\sin (2\pi\,\!x/L),0\right)$
+
+$\mathbf{B}=B_0 \left(-\sin (2\pi\,\!y/L),\sin (4\pi\,\!x/L),0\right)$
+
+$e=p_0/(\gamma -1)+\mathbf{m}^2/(2\varrho)+\mathbf{B}^2/(2\mu)$
+
+where $\varrho_0=2.77778$, $p_0=5/3$ and $B_0=1$. The adiabatic index
+$\gamma =5/3$ and the magnetic permeability $\mu=1$.
 
     void configUser(GRD **gm)
     /**
@@ -1092,22 +1106,30 @@ permeability *μ*=1.
     return;
     }
 
-Example 1 uses the macro `SPQR(a,b,c)` which shortcuts the algebraic
-expression *a*<sup>2</sup> + *b*<sup>2</sup> + *c*<sup>2</sup>.
-
-**Example 2** (taken from `/nirvana/testproblems/MHD/problem17`; cf. \[[Zie05](3.2-User-interfaces#references)\])
-
-IC for the shock-cloud interaction problem simulated in a Cartesian box
-given by (*x*, *y*, *z*) ∈ \[ − 1/2, 1/2\]<sup>3</sup>:
-![ic_example2](uploads/e31ebefeb64b91e0ef74d839ebac4de8/ic_example2.png)  
-At **x**=(0.3,0,0) a spherical clump with radius 0.15 and density of
-10 is embedded and co-moving with its surrounding flow under the
-assumption of pressure equilibrium. The adiabatic index *γ*=5/3 and
-magnetic permeability *μ*=1.
+The macro `SPQR(a,b,c)` defined in header file `nirvana.h` shortcuts the
+algebraic expression $a^2+b^2+c^2$.
+
+**Example 2** (taken from `/nirvana/testproblems/MHD/problem17`; cf.
+[@Z05])
+
+IC for the shock-cloud interaction problem simulated in a 3D Cartesian
+box given by $(x,y,z)\in [-1/2,1/2]^3$:
+$$(\varrho,p,v_x,v_y,v_z,B_x,B_y,B_z)=\left\{\begin{array}{ll}
+(3.86859,167.345,0,0,0,0,2.1826182,-2.1826182) & x<0.1\\
+(1,1,-11.2536,0,0,0,0.56418958,0.56418958) & x\ge 0.1
+\end{array}\right.$$ At ${\mathbf x}=(0.3,0,0)$ a spherical clump with
+radius $0.15$ and density of $10$ is embedded and co-moving with its
+surrounding flow under the assumption of pressure equilibrium. The
+adiabatic index $\gamma =5/3$ and magnetic permeability $\mu=1$. The
+grid is initially refined with 3 refinement levels, `_C.level`=3, in the
+vicinity of the shock and cloud (cf. [User-defined initial mesh
+refinement and refinement
+control](#user-defined-initial-mesh-refinement-and-refinement-control)
+how to refine a grid initially)
 
     void configUser(GRD **gm)
     /**
-     * Definition of ICs.
+     * Definition of IC.
      *
      * @param gm  master mesh pointer
      *
@@ -1205,101 +1227,105 @@ magnetic permeability *μ*=1.
     return;
     }
 
-## Defining boundary conditions
+### Defining boundary conditions
 
 A superblock `g` includes up to 3 ghost cells per side per coordinate
 direction beyond its active region. The number of ghost cells depends on
 the physics and geometry of the problem and is stored in the code
 parameter `_C.ngc`. 3 ghost cells are needed for MHD on
 cylindrical/spherical grids and in simulations with selfgravity, 2 ghost
-cells otherwise.
+cells otherwise. If a superblock shares a physical domain boundary ghost
+cells there have to be assigned with boundary values.
+
+There are 3 categories of boundary conditions (BC): standard,
+user-defined and those for the gravitational potential in simulations
+with selfgravity.
 
 #### Standard BC
 
-There are a number of predefined BC types (standard BC) which can be
-selected in the parameter interface `nirvana.par` under category
-BOUNDARY CONDITION. There, standard types are specified by letters
-{I,O,D,M,A,R,C,F,P} (cf. [Specification of main simulation parameters](3.2-User-interfaces#specification-of-main-simulation-parameters)). 
+Standard BC are predefined and can be selected in the parameter file
+`nirvana.par` below category BOUNDARY CONDITION. Standard types are
+specified by letters out of the set {I,O,D,M,A,R,C,F,P,U}. Each domain
+boundary is assigned a letter. The possible BC types are described in
+detail in [Specification of main simulation
+parameters](#specification-of-main-simulation-parameters). Capital U is
+special and indicates that user-defined BC are to be used at that domain
+boundary.
 
 #### User-defined BC
 
-There is further the option that a user can define its own BC.
-User-defined BC are enabled by specifying letter U in section
-`BOUNDARY CONDITION` in `nirvana.par`. For the actual definition of the
-BC the following user interfaces must be coded, if problem-relevant:
-
-| BC interface  | BC for                                                                     |
-|:--------------|:---------------------------------------------------------------------------|
-| `bcrhoUser.c` | mass density 𝜚                                                             |
-| `bcmUser.c`   | momentum density **m**                                                     |
-| `bceUser.c`   | total energy density *e*                                                   |
-| `bcetUser.c`  | thermal energy density *ε* (only when dual energy mode)                    |
-| `bcbUser.c`   | magnetic field **B**                                                       |
-| `bcnXUser.c`  | species number densities *n*<sub>*s*</sub>, *s* = 0, *N*<sub>*s*</sub> − 1 |
-| `bcCUser.c`   | tracer variables *C*<sub>*c*</sub>, *c* = 0, *N*<sub>*c*</sub> − 1         |
-| `bctfUser.c`  | testfield fluctuations **b**<sub>*t*</sub>, *t* = 0, *N*<sub>*t*</sub> − 1 |
+When U specified BC at the corresponding domain boundary are assumed
+user-defined. In that case the code calls user-programmable functions in
+the following modules:
+
+| module        | function      | BC for                                                                     |
+|:--------------|:--------------|:---------------------------------------------------------------------------|
+| `bcrhoUser.c` | `bcrhoUser()` | mass density 𝜚                                                             |
+| `bcmUser.c`   | `bcmUser()`   | momentum density **m**                                                     |
+| `bceUser.c`   | `bceUser()`   | total energy density *e*                                                   |
+| `bcetUser.c`  | `bcetUser()`  | thermal energy density *ε* (only when dual energy mode)                    |
+| `bcbUser.c`   | `bcbUser()`   | magnetic field **B**                                                       |
+| `bcnXUser.c`  | `bcnXUser()`  | species number densities *n*<sub>*s*</sub>, *s* = 0, *N*<sub>*s*</sub> − 1 |
+| `bcCUser.c`   | `bcCUser()`   | tracer variables *C*<sub>*c*</sub>, *c* = 0, *N*<sub>*c*</sub> − 1         |
+| `bctfUser.c`  | `bctfUser()`  | testfield fluctuations **b**<sub>*t*</sub>, *t* = 0, *N*<sub>*t*</sub> − 1 |
 
 A combination of user-defined BC and standard BC at different domain
-boundaries is possible. Examples can be found in the testproblems
-collection located in subdirectories under `/nirvana/testproblems/`,
+boundaries is possible. Examples can be found in NIRVANA's collection of
+testproblems located in subdirectories below `/nirvana/testproblems/`,
 e.g., the double Mach reflection problem defined in
-`/nirvana/testproblems/MHD/problem1`.
+`/nirvana/testproblems/MHD/problem1`. Note that there is no need to
+define BC for variables which are not used by the problem.
+
+**Important:** Although not a primary variable, defining BC for the
+thermal energy density in `bcetUser()` is mandatory if the dual energy
+formalism is enabled (cf. [Specification of main simulation
+parameters](#specification-of-main-simulation-parameters)).
 
-At the *x*(*y*,*z*)-lower domain boundary of a superblock `g` boundary
-values have to be assigned from `ix`=0 to `g->ixs`-1 (`iy`=0 to
-`g->iys`-1, `iz`=0 to `g->izs`-1) and full range for the cross indices.
-Recall that `g->ixs` (`g->iys`, `g->izs`) counts the first active cell
-on `g` in *x*(*y*,*z*)-direction. At the *x*(*y*,*z*)-upper domain
-boundary the `ix`(`iy`,`iz`)-index range of ghost cells depends on the
-physical variable as per
+Boundary values have to be assigned in ghost cells running from to
+`g->ixs`-1 ( to `g->iys`-1, to `g->izs`-1) in $x$($y$,$z$)-direction at
+the lower domain boundary. At upper domain boundaries the index range
+depends on the type of variable as listed in the following table:
 
-| variable                                                                                                   | ix-range (x-upper)   | iy-range (y-upper)   | iz-range (z-upper)        |
+| variable                                                                                                   | ix-range             | iy-range             | iz-range                  |
 |:-----------------------------------------------------------------------------------------------------------|:---------------------|:---------------------|:--------------------------|
 | 𝜚, *m*<sub>*x*</sub>, *m*<sub>*y*</sub>, *m*<sub>*z*</sub>, *e*, *ε*, *n*<sub>*s*</sub>, *C*<sub>*c*</sub> | `g->ixe`+1,`g->nx`   | `g->iye`+1,`g->ny`   | `g->ize`+`_C.idz`,`g->nz` |
 | *B*<sub>*x*</sub>,*b*<sub>*x*<sub>*t*</sub></sub>                                                          | `g->ixe`+2,`g->nx`+1 | `g->iye`+1,`g->ny`   | `g->ize`+`_C.idz`,`g->nz` |
 | *B*<sub>*y*</sub>,*b*<sub>*y*<sub>*t*</sub></sub>                                                          | `g->ixe`+1,`g->nx`   | `g->iye`+2,`g->ny`+1 | `g->ize`+`_C.idz`,`g->nz` |
 | *B*<sub>*z*</sub>,*b*<sub>*z*<sub>*t*</sub></sub>                                                          | `g->ixe`+1,`g->nx`   | `g->iye`+1,`g->ny`   | `g->ize`+2,`g->nz`+1      |
 
-and cross indices again run through their full range. For instance, the
-following code fragment is the loop for assigning the *x*-upper BC for
-*B*<sub>*x*</sub>:
-
-      /* UPPER X-FACE */
-      for(ix=g->ixe+2; ix<=g->nx+1; ix++) 
-        for(iz=0; iz<=g->nz; iz++) 
-          for(iy=0; iy<=g->ny; iy++) 
-          {
-            g->bx[iz][iy][ix]=...
-          }
-
 #### BC for the gravitational potential
 
-BC for the gravitational potential *Φ* cannot be explicitely set.
-Instead, the following rule is adopted which transforms speficied types
-into BC types for *Φ*:
+BC for the gravitational potential $\Phi$ cannot be explicitely set
+except when specifying U in `nirvana.par`. Otherwise, the following rule
+is adopted to transform MHD BC types into BC types for $\Phi$:
 
 | selected BC type |  ⇒  | BC type for *Φ*                        |
-|-----------------:|:---:|:----------------|
-|                P |     | P               |
+|------------------|:----|:---------------------------------------|
+|                P |     | P (periodicity)                        |
 |            M,A,R |     | von-Neumann                            |
 |            I,O,D |     | Dirichlet                              |
-| user-defined (U) |     | Dirichlet       |
+|                F |     | not supported                          |
+|                U |     | user-defined *Φ* (function `phiUser()` |
+
+In case of von-Neumann conditions (M,A,R) the gradient of potential
+vanishes normal to the corresponding domain boundary, i.e.,
+$\mathbf{e}_n\cdot \nabla\Phi=0$ where $\mathbf{e}_n$ is the unit normal
+vector of the domain boundary. In case of Dirichlet conditions (I,O,D)
+boundary values for $\Phi$ are obtained by a multipole expansion for the
+given mass distribution (cf. [physics
+guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
+In case of user-defined conditions U the selfgravity solver calls the
+function `phiUser()` expecting a user-defined gravitational potential.
 
-In case of von-Neumann conditions the gradient of potential vanishes
-normal to the corresponding domain boundary, i.e.,
-**e**<sub>*n*</sub> ⋅ ∇*Φ* = 0 where **e**<sub>*n*</sub> is the unit
-normal vector of the domain boundary. In case of Dirichlet conditions
-boundary values for *Φ* are obtained by a multipole expansion for the
-given mass distribution. More details and formulae can be found in the
-physics guide
-([Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/A-PhysicsGuide),
-[PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/PhysicsGuide.pdf)).
+**Important:** A combination of periodic BC with non-periodic BC is not
+possible. The solver only supports triple-periodic BC when periodic.
 
-*Note: The selfgravity solver supports only triple-periodic BC or
-non-periodic BC. The combination of periodic BC at one domain boundary
-with non-periodic BC at another domain boundary is not supported.*
+When using the screening mass approach for cylindrical/spherical
+problems with $2\pi$-periodicity BC are computed to discretization error
+accuracy by the gravity solver itself. In this case above rules do not
+apply and any specified BC are ignored.
 
-## User-defined coefficients for dissipative processes
+### User-defined coefficients for dissipative processes
 
 The modules `viscosityCoeffUser.c`, `conductionCoeffUser.c` and
 `diffusionCoeffUser.c` serve as templates for a user-defined coefficient
@@ -1309,28 +1335,30 @@ respectively.
 #### Fluid viscosity
 
 Fluid viscosity enters the momentum equation and energy equation via the
-viscous stress tensor *τ* given by
+viscous stress tensor $\tau$ given by
 
-*τ* = *ν*\[∇**v**+(∇**v**)<sup>⊤</sup>−2/3(∇⋅**v**)*I*\]
+$$\tau=\nu \left[\nabla\mathbf{v}+(\nabla\mathbf{v})^\top
+-2/3(\nabla\cdot\mathbf{v})\mathrm{I}\right]$$
 
-where *ν* \[kg⋅m<sup>−1</sup>⋅s<sup>−1</sup>\] is the
-*dynamic* viscosity coefficient, **v** the fluid velocity and *I* the
-identity operator.
+where $\nu$ in units
+\[$\mathtt{kg}\cdot\mathtt{m}^{-1}\cdot\mathtt{s}^{-1}$\] is the
+*dynamic* viscosity coefficient, $\mathbf{v}$ the fluid velocity and
+$\mathrm{I}$ the identity operator.
 
-*Note: The *dynamic* viscosity coefficient *ν* to be defined here should
-not be confused with the kinetic coefficient,*
-*ν*<sub>*kinetic*</sub>
-\[m<sup>2</sup>⋅s<sup>−1</sup>\], *related to the dynamic
-coefficient by* *ν* = 𝜚*ν*<sub>*kinetic*</sub>.
+**Important:** The viscosity coefficient $\nu$ here should not be
+confused with the kinetic viscosity coefficient,
+$\nu_{\mathrm{kinetic}}$ with unit \[m$^2\cdot$s$^{-1}$\], which is
+related to the dynamic coefficient by
+$\nu=\varrho\nu_{\mathrm{kinetic}}$.
 
-In the call of `viscosityCoeffUser()` the function arguments are the
-superblock pointer `g` and the array pointer `vis` of type `double***`
-representing the viscosity coefficient *ν*:
+A user-defined coefficient, $\nu$, has to be assigned in module
+`viscosityUser.c` in the function
 
       viscosityCoeffUser(g,vis);
 
-The user must assign the `vis`-array for the full index range of `g`,
-whence a loop like
+taking arguments `g`, a superblock pointer (not the mesh master
+pointer), and the array pointer `vis` of type `double***` representing
+the dynamic viscosity coefficient. Here is an example code fragment
 
       /* MY VISCOSITY COEFFICIENT */
       for(iz=0; iz<=g->nz; iz++)
@@ -1340,41 +1368,42 @@ whence a loop like
             vis[iz][iy][ix]=...
           }
 
-must be added in `viscosityCoeffUser.c`. The coefficient should be
-evaluated at cell-centroid coordinates
-(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). An example definition for a viscosity
-coefficient can be found in testproblem
-`/nirvana/testproblems/VISC/problem1`.
+Spatial indexing covers here the full domain of the superblock including
+ghost cells. The coefficient is assumed a cell-centroid quantity, i.e.,
+it should be defined at location (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`).
 
-User-defined viscosity is enabled by appropriate choice in the parameter
-interface `nirvana.par` under the category PHYSICS SPECIFICATIONS
-(code parameter: `_C.viscosity`).
+Fluid viscosity with user-defined viscosity coefficient is enabled by
+the appropriate choice of parameter `_C.viscosity` in file `nirvana.par`
+below category PHYSICS SPECIFICATIONS. An example can be found in
+testproblem `/nirvana/testproblems/VISC/problem1`.
 
 #### Thermal conduction
 
 Thermal conduction enters the energy equation through a heat flux
-**F**<sub>*C*</sub>. Generally, the presence of a magnetic field
+$\mathbf{F}_{\mathrm{C}}$. Generally, the presence of a magnetic field
 introduces anisotropic effects with different transport properties along
-and across the magnetic field. **F**<sub>*C*</sub> is described by
+and across the magnetic field. $\mathbf{F}_\mathrm{C}$ is described by
 
-**F**<sub>*C*</sub> = −*κ*<sub>∥</sub>(∇*T*⋅**B̂**)**B̂**−*κ*<sub>⊥</sub>(∇*T*−(∇*T*⋅**B̂**)**B̂**)
+$$\mathbf{F}_\mathrm{C}=-\kappa_\parallel (\nabla T\cdot\mathbf{\hat{B}})\mathbf{\hat{B}}
+-\kappa_\perp \left(\nabla T-(\nabla T\cdot\mathbf{\hat{B}})\mathbf{\hat{B}}\right)$$
 
-where *κ*<sub>∥</sub> and *κ*<sub>⊥</sub> is the thermal conduction
+where $\kappa_\parallel$ and $\kappa_\perp$ is the thermal conduction
 coefficient parallel and perpendicular to the magnetic field,
-respectively, and **B̂** = **B**/\|**B**\| is the unit vector in the
-direction of the magnetic field. *κ*<sub>∥</sub> and *κ*<sub>⊥</sub> are
-measured in units
-J⋅K<sup>−1</sup>⋅m<sup>−1</sup>⋅s<sup>−1</sup>.
+respectively, and $\mathbf{\hat{B}}=\mathbf{B}/|\mathbf{B}|$ is the unit
+vector in the direction of the magnetic field. $\kappa_\parallel$ and
+$\kappa_\perp$ are measured in units
+J$\cdot$K$^{-1}\cdot$m$^{-1}\cdot$s$^{-1}$ (cf. [physics
+guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
 
-In the call of `conductionCoeffUser()` the function arguments are the
-superblock pointer `g` and the array pointers `cond`, `cond_perp` of
-type `double***` representing the conduction coefficients
-*κ*<sub>∥</sub>, *κ*<sub>⊥</sub>:
+User-defined coefficients, $\kappa_\parallel$ and $\kappa_\perp$, have
+to be assigned in module `conductionCoeffUser.c` in the function
 
       conductionCoeffUser(g,cond,cond_perp);
 
-The user must assign the `cond,cond_perp`-arrays for the full index
-range of `g`, whence a loop like
+taking arguments `g`, a superblock pointer (not the mesh master
+pointer), and the array pointers `cond` and `cond_perp` of type
+`double***` representing the parallel and perpendicular conduction
+coefficients, respectively. Here is an example code fragment:
 
       /* MY CONDUCTION COEFFICIENT */
       for(iz=0; iz<=g->nz; iz++)
@@ -1386,36 +1415,37 @@ range of `g`, whence a loop like
             cond_perp[iz][iy][ix]=...
           }
 
-must be added to `conductionCoeffUser.c`. The coefficients should be
-evaluated at cell-centroid coordinates
-(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). In HD simulations or MHD simulations
-with forced isotropy (when `COND_FORCE_ISO`=`YES` in `nirvanaUser.h`)
-only the `cond`-array has to be assigned with `cond` representing the
-conduction coefficient *κ* in the isotropic heat flux
-**F**<sub>*C*</sub> = −*κ*∇*T*.
+Spatial indexing covers here the full domain of the superblock including
+ghost cells. Both coefficients are assumed cell-centroid quantities,
+i.e., it should be defined at location
+(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). In HD simulations and in MHD
+simulations with forced isotropy (when `COND_FORCE_ISO`=`YES` in
+`nirvanaUser.h`) only the `cond`-array has to be assigned with `cond`
+representing now the conduction coefficient $\kappa$ in the isotropic
+heat flux $\mathbf{F}_\mathrm{C}=-\kappa \nabla T$.
 
-User-defined thermal conduction is enabled by appropriate choice in the
-parameter interface `nirvana.par` under the category
-PHYSICS SPECIFICATIONS (code parameter: `_C.conduction`).
+Thermal conduction with user-defined conduction coefficients is enabled
+by the appropriate choice of parameter `_C.conduction` in file
+`nirvana.par` below category PHYSICS SPECIFICATIONS.
 
 #### Ohmic diffusion
 
 Ohmic diffusion enters the induction equation and energy equation as a
 field contribution given by
 
-**E**<sub>*D*</sub> = *η*<sub>*D*</sub>∇×**B**
+$\mathbf{E}_\mathrm{D}=\eta_\mathrm{D} \nabla\times\mathbf{B}$
 
-where *η*<sub>*D*</sub> \[m<sup>2</sup>⋅s<sup>−1</sup>\] is the
-diffusion coefficient.
+where $\eta_\mathrm{D}$ in units \[m$^2\cdot$s$^{-1}$\] is the diffusion
+coefficient.
 
-In the call of `diffusionCoeffUser()` function arguments are the
-superblock pointer `g` and the array pointer `diff` of type `double***`
-representing the diffusion coefficient *η*<sub>*D*</sub>:
+A user-defined coefficients, $\eta_\mathrm{D}$, has to be assigned in
+module `diffusionCoeffUser.c` in the function
 
       diffusionCoeffUser(g,diff);
 
-The user must assign the `diff`-array for the full index range of `g`,
-whence a loop like
+taking arguments `g`, a superblock pointer (not the mesh master
+pointer), and the array pointer `diff` of type `double***` representing
+the diffusion coefficient. Here is an example code fragment:
 
       /* MY DIFFUSION COEFFICIENT */
       for(iz=0; iz<=g->nz; iz++)
@@ -1425,39 +1455,39 @@ whence a loop like
             diff[iz][iy][ix]=...
           }
 
-must be added in `diffusionCoeffUser.c`. The coefficient should be
-evaluated at cell-centroid coordinates
+Spatial indexing covers here the full domain of the superblock including
+ghost cells. The diffusion coefficient is assumed a cell-centroid
+quantity, i.e., it should be defined at location
 (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`).
 
-User-defined Ohmic diffusion is enabled by appropriate choice in the
-parameter interface `nirvana.par` under the category
-PHYSICS SPECIFICATIONS (code parameter: `_C.diffusion`).
+Ohmic diffusion with user-defined diffusion coefficient is enabled by
+appropriate choice of the parameter `_C.diffusion` in `nirvana.par`
+below category PHYSICS SPECIFICATIONS.
 
-## User-defined coefficient for ambipolar diffusion
+### User-defined coefficient for ambipolar diffusion
 
 The module `APdiffusionCoeffUser.c` serves as template for a
-user-defined ambipolar diffusion coefficient. Ambipolar diffusion enters
-the induction equation and energy equation as a field contribution given
-by
+user-defined ambipolar diffusion coefficient.
 
-**E**<sub>*AD*</sub> = *η*<sub>*AD*</sub>/*μ* **B**×[(∇×**B**)×**B**\]
+Ambipolar diffusion enters the induction equation and energy equation as
+a field contribution given by
 
-where *η*<sub>*AD*</sub>
-\[V⋅m⋅A<sup>−1</sup>⋅T<sup>−2</sup>\] denotes the
-ambipolar diffusion coefficient.
+$\mathbf{E}_\mathrm{AD}=\eta_\mathrm{AD}/\mu\mathbf{B}\times\left[(\nabla\times\mathbf{B})
+\times\mathbf{B}\right]$
 
-*Note: The prefactor* *η*<sub>*AD*</sub>/*μ* *has units*
-m<sup>2</sup>⋅s<sup>−1</sup>⋅T<sup>−2</sup>.
+where $\eta_\mathrm{AD}$ \[V$\cdot$m$\cdot$A$^{-1}\cdot$T$^{-2}$\]
+denotes the ambipolar diffusion coefficient. The prefactor
+$\eta_\mathrm{AD}/\mu$ has units m$^2\cdot$s$^{-1}\cdot$T$^{-2}$.
 
-In the call of `APdiffusionCoeffUser()` function arguments are the
-superblock pointer `g` and the array pointer `APdiff` of type
-`double***` representing the ambipolar diffusion coefficient
-*η*<sub>*AD*</sub>:
+A user-defined coefficients, $\eta_\mathrm{AD}$, has to be assigned in
+module `APdiffusionCoeffUser.c` in the function
 
       APdiffusionCoeffUser(g,APdiff);
 
-The user must assign the `APdiff`-array for the full index range of `g`,
-whence a loop like
+taking arguments `g`, a superblock pointer (not the mesh master
+pointer), and the array pointer `APdiff` of type `double***`
+representing the ambipolar diffusion coefficient. Here is an example
+code fragment:
 
       /* MY AMBIPOLAR DIFFUSION COEFFICIENT */
       for(iz=0; iz<=nz; iz++)
@@ -1467,31 +1497,30 @@ whence a loop like
             APdiff[iz][iy][ix]=...
           }
 
-must be added in `APdiffusionCoeffUser.c`. The coefficient should be
-evaluated at cell-centroid coordinates
-(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). An example definition for an
-ambipolar diffusion coefficient can be found in testproblem
-`/nirvana/testproblems/APDIFF/problem1`.
+Spatial indexing covers here the full domain of the superblock including
+ghost cells. The ambipolar diffusion coefficient is assumed a
+cell-centroid quantity, i.e., it should be defined at location
+(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`).
 
-User-defined ambipolar diffusion is enabled by appropriate choice in the
-parameter interface `nirvana.par` under the category
-PHYSICS SPECIFICATIONS (code parameter: `_C.APdiffusion`).
+Ambipolar diffusion with user-defined ambipolar diffusion coefficient is
+enabled by appropriate choice of the parameter `_C.APdiffusion` in
+`nirvana.par` below category PHYSICS SPECIFICATIONS. An example can be
+found in testproblem `/nirvana/testproblems/APDIFF/problem1`.
 
-## User-defined body force
+### User-defined body force
 
-The module `forceUser.c` serves as template for coding a user-defined
-*specific* body force **f** (force per mass in units
-N⋅kg<sup>−1</sup>). The body force then enters the momentum
-equation and energy equation as source term.
+The module `forceUser.c` serves as template for a user-defined
+*specific* body force $\mathbf{f}$ (force per mass in units
+N$\cdot$kg$^{-1}$). The body force enters the momentum equation and
+energy equation as source term.
 
-In the call of `forceUser()` the function arguments are the superblock
-pointer `g` and the array pointers `fx`,`fy`,`fz` of type `double***`
-representing the force components:
+A user body force has to be defined in the function
 
       forceUser(g,fx,fy,fz);
 
-The user must assign the `fx`,`fy`,`fz`-arrays for active grid cells of
-`g`, whence a loop like
+taking arguments `g`, a superblock pointer (not the mesh master
+pointer), and the array pointers `fx`,`fy`,`fz` of type `double***`
+representing the force components. Here is an example code fragment:
 
       /* MY BODY FORCE */
       for(iz=g->izs; iz<=g->ize; iz++)
@@ -1505,76 +1534,62 @@ The user must assign the `fx`,`fy`,`fz`-arrays for active grid cells of
             fz[iz][iy][ix]=...
           }
 
-must be added in `forceUser.c`. Strictly speaking, the body force is to
-be understood as a cell-averaged quantity. Therefore, **f** should be
-evaluated at cell-centroid coordinates (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`)
-and the point values assigned to `fx`,`fy`,`fz`. An example definition
-for a body force can be found in testproblem
-`/nirvana/testproblems/MHD/problem21`.
+Spatial indexing covers the active domain of the superblock. The force
+components are assumed cell-centroid quantities, i.e., it should be
+defined at location (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`).
 
-The body force is enabled by specification in the parameter interface
-`nirvana.par` under the category PHYSICS SPECIFICATIONS (code
-parameter: `_C.force`).
+A body force is enabled by specification of the parameter `_C.force` in
+file `nirvana.par` below category PHYSICS SPECIFICATIONS. An example can
+be found in testproblem `/nirvana/testproblems/MHD/problem21`.
 
-## User-defined cooling and heating function
+### User-defined cooling- and heating function
 
 The modules `sourceCoolingUser.c` and `sourceHeatingUser.c` serve as
 templates for coding a user-defined cooling function,
-*L*<sub>*cool*</sub>≤0, and heating function,
-*L*<sub>*heat*</sub>≥0, respectively. Both functions are allowed
-to depend on temperature *T* and density 𝜚. The net heatloss, i.e. the
-sum
-*L*<sub>*cool*</sub>(*T*,𝜚)+*L*<sub>*heat*</sub>(*T*,𝜚)
-of both functions, enters as a source term in the energy equation.
-*L*<sub>*cool*</sub> and *L*<sub>*heat*</sub> are measured
-in units J⋅s<sup>−1</sup>⋅m<sup>−3</sup>.
-
-In the call of `sourceCoolingUser()` (`sourceHeatingUser()`) the
-function arguments are the temperature value `T`, density value `rho`,
-the pointer `deriv` to the derivatives flag and the 2-element vector
-`dfc` (`dfh`):
-
-      f=sourceCoolingUser(T,rho,deriv,dfc);
-
-      f=sourceHeatingUser(T,rho,deriv,dfh);
-
-The user must define the return value
-`f` = *L*<sub>*cool*</sub>(`T`,`rho`)
-(*L*<sub>*heat*</sub>(`T`,`rho`)) in `sourceCoolingUser.c`
-(`sourceHeatingUser.c`). The `deriv`-flag is thought to indicate the
-calling function whether a user provides the derivatives of
-*L*<sub>*cool*</sub>(*T*,𝜚) and
-*L*<sub>*heat*</sub>(*T*,𝜚) with respect to *T* and 𝜚 himself.
-The `deriv`-flag is preset to `NO` when entering the code function. If
-the user sets
-
-      *deriv=YES; 
-
-the derivatives must be stored by the user in the 2-element vector `dfc`
-(`dfh`) with `dfc[0]` (`dfh[0]`) storing the derivative
-∂*L*<sub>*cool*</sub>/∂*T* (∂*L*<sub>*heat*</sub>/∂*T*) and
-`dfc[1]` (`dfh[1]`) storing the derivative
-∂*L*<sub>*cool*</sub>/∂𝜚 (∂*L*<sub>*heat*</sub>/∂𝜚).
-Otherwise derivatives are automatically computed by a finite difference
-approximation.
-
-*Note: Knowledge of the derivatives is a prerequisite in the numerical
-procedure, the Field-length-based mesh refinement criterion and timestep
-computation.*
-
-An example definition for a cooling/heating function can be found in
-testproblem `/nirvana/testproblems/HEATLOSS/problem2`.
-
-User-defined cooling/heating is enabled by appropriate choice in the
-parameter interface `nirvana.par` under the category
-PHYSICS SPECIFICATIONS (code parameter: `_C.heatloss`).
-
-## User-defined equation of state
+$L_\mathrm{cool}\le 0$, and heating function, $L_\mathrm{heat}\ge 0$,
+respectively. Both functions are allowed to depend on temperature $T$
+and density $\varrho$. The net heatloss, i.e. the sum
+$L_\mathrm{cool}(T,\varrho)+L_\mathrm{heat}(T,\varrho)$, enters as
+source term in the energy equation. $L_\mathrm{cool}$ and
+$L_\mathrm{heat}$ are measured in units J$\cdot$s$^{-1}\cdot$m$^{-3}$.
+
+User-defined cooling/heating rates are to be defined in the functions
+
+      f=sourceCoolingUser(T,rho,deriv,dfc)
+
+      f=sourceHeatingUser(T,rho,deriv,dfh)
+
+taking arguments `T`, the temperature, `rho`, the gas density, the
+pointer `deriv`, a flag to tell the calling function whether the user
+provides the derivatives $\partial L_\mathrm{cool}\partial\varrho$,
+$\partial L_\mathrm{cool}\partial T$ to be stored in the 2-element
+vector `dfc[0], dfc[1]` and the derivatives
+$\partial L_\mathrm{heat}\partial\varrho$,
+$\partial L_\mathrm{heat}\partial T$ to be stored in the 2-element
+vector `dfh[0], dfh[1]`. The default value is preset to `deriv=YES` in
+both functions. The user may avoid to explicitely calculating
+derivatives by setting `deriv=NO`. Then, derivatives are computed by the
+code with a finite difference approximation. The functions return `f`,
+the user-defined cooling/heating rates.
+
+**Important:** Defining derivatives of the cooling/heating rates is a
+prerequisite for the exponential Rosenbrock numerical method, the
+Field-length-based mesh refinement criterion and for the cooling/heating
+timestep computation.
+
+User-defined cooling/heating is enabled by appropriately choosing the
+parameter `_C.heatloss` in file `nirvana.par` below category PHYSICS
+SPECIFICATIONS. An example can be found in testproblem
+`/nirvana/testproblems/HEATLOSS/problem2`.
+
+### User-defined equation of state
+
+There are two types of an user-defined EOS: analytic and tabulated.
 
 #### Analytic EOS
 
-NIRVANA makes possible a simple user-defined implementation of an
-analytic EOS. This requires the defintion of the macros
+NIRVANA allows a simple user-defined implementation of an analytic EOS
+by definition of the following macros:
 
 | macro             | thermodynamic quantity |
 |:------------------|:-----------------------|
@@ -1583,15 +1598,16 @@ analytic EOS. This requires the defintion of the macros
 | `TUSR(rho,eth)`   | temperature            |
 | `ETUSR(rho,eth)`  | thermal energy density |
 
-as a function of mass density 𝜚 (`rho`) and thermal energy density *ε*
-(`eth`). The macros are declared in the header file `nirvanaUser.h`. The
-following code fragment (taken from testproblem
-`/nirvana/testproblems/GRAVITY/problem4`) illustrates the definition of
-a barotropic EOS of the form
+The macros are in general functions of the mass density $\varrho$
+(`rho`) and thermal energy density $\varepsilon$ (`eth`). It are
+declared in the header file `nirvanaUser.h`. In order to give an
+example, the following code fragment (taken from `nirvanaUser.h` in
+testproblem `/nirvana/testproblems/GRAVITY/problem4`) illustrates the
+definition of a barotropic EOS of the form
 
-*p*(𝜚) = *a*<sup>2</sup>𝜚\[1+(𝜚/𝜚<sub>0</sub>)<sup>4/3</sup>\]<sup>1/2</sup>
+$p(\varrho)=a^2\varrho \left[1+(\varrho/\varrho_0)^{4/3}\right]^{1/2}$
 
-with *a*=200 and 𝜚<sub>0</sub>=10<sup>−10</sup>.
+with $a=200$ and $\varrho_0=10^{-10}$.
 
     #define PUSR(rho,eth)   (40000.*(rho)*sqrt(1.+pow((rho)/1.e-10,4./3.))+0*(eth))
     #define CS2USR(rho,eth) (40000.*sqrt(1.+pow((rho)/1.e-10,4./3.)) \
@@ -1600,78 +1616,82 @@ with *a*=200 and 𝜚<sub>0</sub>=10<sup>−10</sup>.
     #define TUSR(rho,eth)   (MOK*_C.mean_molecular_weight*PUSR(rho,eth)/rho)
     #define ETUSR(rho,eth)  (PUSR(rho,eth)/(_C.gamma-1.))
 
-Since *p* solely depends on 𝜚 the energy density parameter `eth` is
-irrelevant is this case. The `ETUSR` definition makes use of the prior
-definition of `PUSR` and assumes the usual caloric EOS with ratio of
-specific heats *γ* (code variable: `_C.gamma`). No energy equation is
-solved in this testproblem.
+Since $p$ solely depends on $\varrho$ the energy density parameter `eth`
+is irrelevant is this case. The `ETUSR` definition makes use of the
+prior definition of `PUSR` and assumes the usual caloric EOS with ratio
+of specific heats $\gamma$ (code variable `_C.gamma`). A barotropic EOS
+is consistent with solving no energy equation.
 
-*Note 1: If the energy equation is solved together with an analytic EOS
-the macro `ETUSR` must read*
+**Important:** If an energy equation is solved together with an analytic
+EOS the macro `ETUSR` must necessarily read
 
     #define ETUSR(rho,eth)  (eth)
 
-*for obvious reasons of consistency.*
+for obvious reasons of consistency.
 
-*Note 2: The macros are used in function `TDeos()` in module `utilTD.c`
-and are needed in the numerical procedure.*
+**Important:** In case of an analytic EOS the macros are used in the
+function `TDeos()` in module `utilTD.c`, and it are also used at other
+places in the numerical solver.
 
-An analytic EOS is enabled by appropriate choice in the parameter
-interface `nirvana.par` under category PHYSICS SPECIFICATIONS (code
-paramter: `_C.eos`).
+An analytic EOS is enabled by appropriately choosing the parameter
+`_C.eos` in file `nirvana.par` below category PHYSICS SPECIFICATIONS.
 
 #### Tabulated EOS
 
-NIRVANA also permits the implementation of a tabulated EOS. It requires
-the creation of look-up tables for the logarithm of pressure, log *p*,
-and logarithm of temperature, log *T*, both as functions of the
-logarithm of mass density 𝜚 and the logarithm of thermal energy density
-*ε*. The look-up tables for a given (log 𝜚, log *ε*)-domain are
-generated by cubic interpolation from user-specified sampling data. The
-sampling data {(log *p*)<sub>*ij*</sub>}
-({(log *T*)<sub>*ij*</sub>}) on the grid
-{(log 𝜚)<sub>*i*</sub>} × {(log *ε*)<sub>*j*</sub>} has to be defined by
-the user in the interface `eosTabPressureUser.c`
-(`eosTabTemperatureUser.c`). The discretizations in log 𝜚-space and
-log *ε*-space, {(log 𝜚)<sub>*i*</sub>; i=0,`n1`-1} and {(log *ε*)<sub>*j*</sub>; j=0,`n2`-1}, must be equidistant
-where `n1` and `n2` are the number of sampling points.
-
-Actually, in `eosTabPressureUser()` (`eosTabTemperatureUser()`) the user
-must return a structure of type `UTAB` (declared in header `tabular.h`)
-which collects all relevant sampling data. Denoting by `ut` an instance
-of `UTAB` the structure elements are accessed by `ut.V`. The following
-structure elements `V` must be defined by the user:
+NIRVANA also permits the implementation of a tabulated EOS. However, the
+procedure is more involved and requires the creation of look-up tables
+for the logarithm of pressure, $\log p$, and logarithm of temperature,
+$\log T$, both as functions of the logarithm of mass density $\varrho$
+and the logarithm of thermal energy density $\varepsilon$. The look-up
+tables for a prescribed $(\log\varrho, \log\varepsilon)$-range are
+generated by cubic interpolation from user-specified data
+$\{(\log p)_{ij}\}$ and $\{(\log T)_{ij}\}$ on an rectilinear mesh
+$\{(\log\varrho)_i\}\times\{(\log\varepsilon)_j\}$, $i=0,{\tt n1}-1$ and
+$j=0,{\tt n2}-1$ with `n1`$\times$`n2` the number of mesh points.
+
+Data for the look-up tables has to be defined in the function
+`eosTabPressureUser()` in module `eosTabPressureUser.c` for $\log p$,
+and in the function `eosTabTemperatureUser()` in module
+`eosTabTemperatureUser.c` for $\log T$.
+
+In both, `eosTabPressureUser()` and `eosTabTemperatureUser()`, the user
+must assign and return a structure of type `UTAB` (declared in header
+`tabular.h`) which collects all relevant data to construct the tables.
+Denoting by `ut` an instance of `UTAB` the structure is accessed by
+`ut.V` with elements `V` summarized in the following table:
 
 | `V`             | meaning                                                                                |
 |:----------------|:---------------------------------------------------------------------------------------|
-| `n1`            | number of sampling points in log 𝜚                                                     |
-| `n2`            | number of sampling points in log *ε*                                                   |
+| `n1`            | number of mesh points in log 𝜚                                                     |
+| `n2`            | number of mesh points in log *ε*                                                   |
 | `u1_lo`,`u1_up` | lower,upper value of log 𝜚-range                                                       |
 | `u2_lo`,`u2_up` | lower,upper value of log *ε*-range                                                     |
-| `dt2[i][j]`     | 2D array for sampling data {(log *p*)<sub>*ij*</sub>} ({(log *T*)<sub>*ij*</sub>}) |
+| `dt2[i][j]`     | 2D data array for {(log *p*)<sub>*ij*</sub>} ({(log *T*)<sub>*ij*</sub>}) |
 
-The subsequent example taken from testproblem
-`/nirvana/testproblem/MHD/problem1B` provides sampling data for
-tabulation of the caloric EOS *p*=(*γ* − 1)*ε*.
+Here is an example how to code `eosTabPressureUser()` taken from
+testproblem `/nirvana/testproblem/MHD/problem1B` which provides data for
+the standard EOS of caloric type, $p=(\gamma -1)\varepsilon$, in
+tabulated form:
 
     UTAB eosTabPressureUser(void)
     /**
-     * Template for the definition of user-defined sampling data for ln(pressure)
-     * on 2D [ln(density),ln(thermal_energy_density)]-space. The sample data is 
-     * used for creation of the 2D look-up table _TABLP in a tabulated EOS.
+     * Template for the definition of user-defined data for ln(pressure) as a
+     * function of (ln(density),ln(thermal_energy_density)). The returned structure 
+     * is used to create the look-up table _TABLP.
      *
-     * @retval ut  UTAB structure for creating the 2D look-up table _TABLP
+     * @retval ut  UTAB structure for creating the look-up table _TABLP
      *
      * @author  Udo Ziegler, AIP
      *
      * @note  The UTAB structure is declared in tabular.h and contains the 
      *        following elements which has to be assigned here by the user:
-     *                        ut.n1 = number of equidistant points for ln(density);
+     *                        ut.n1 = number of equidistant points for ln(density)
      *                        ut.n2 = number of equidistant points for 
-     *                                ln(thermal_energy_density);
-     *          [ut.u1_lo,ut.u1_up] = ln(density) variable range;
-     *          [ut.u2_lo,ut.u2_up] = ln(thermal_energy_density) variable range;
-     *                       ut.dt2 = 2D array for ln(pressure) sampling data.
+     *                                ln(thermal_energy_density)
+     *          [ut.u1_lo,ut.u1_up] = ln(density) range
+     *          [ut.u2_lo,ut.u2_up] = ln(thermal_energy_density) range
+     *                       ut.dt2 = 2D data array for ln(pressure)
+     *
      *        Example implementation: tabulated ideal EOS for testproblem 
      *                                MHD/problem1B.
      */
@@ -1702,7 +1722,7 @@ tabulation of the caloric EOS *p*=(*γ* − 1)*ε*.
     for(i2=1; i2<n2; i2++) 
       let[i2]=let[i2-1]+dlet;
 
-    /* LN(PRESSURE) SAMPLING DATA FOR ln(p)=ln(_C.gamma-1)+ln(et) */
+    /* LN(PRESSURE) DATA FOR ln(p)=ln(_C.gamma-1)+ln(et) */
     ut.dt2=Array2(n1-1,n2-1);
 
     lc=log(_C.gamma-1.);
@@ -1715,57 +1735,56 @@ tabulation of the caloric EOS *p*=(*γ* − 1)*ε*.
     return ut;
     }
 
-First, the number of sampling points `ut.n1` in (log 𝜚)-direction and
-`ut.n2` in (log *ε*)-direction are set followed by the ranges
-\[`ut.u1_lo`,`ut.u1_up`\] for log 𝜚 and \[`ut.u2_lo`,`ut.u2_up`\] for
-log *ε*, respectively. Then, the 2D array `ut.dt2` is allocated and
-sampling data {(log *p*)<sub>*ij*</sub>} is assigned. The returned
-`ut` is used by the calling function to finally create the public
-look-up table `_TABLP`. A corresponding user definition for log *T* in
-`eosTabTemperatureUser.c` provides the input for the public look-up
-table `_TABLT`.
-
-*Note: `_TABLP` and `_TABLT` are actually pointers of struct type `TAB`
-(declared in the header `tabular.h`). It are used in the function
-`TDeos()` in module `utilTD.c` and in the numerical procedure.*
-
-A tabulated EOS is enabled by appropriate choice in the parameter
-interface `nirvana.par` under category PHYSICS SPECIFICATIONS (code
-paramter: `_C.eos`).
-
-## User-defined mesh refinement
-
-For certain problems it may be advantageous to start a simulation with a
-pre-refined mesh in some parts of the computational domain or to retrict
-mesh refinement to some parts of it during runtime. Such situations can
-be realize with help of the user interfaces `initDomainUser.c` and
+First, the number of mesh points, `ut.n1` and `ut.n2`, in
+$\log\varrho$-direction and $\log\varepsilon$-direction is set. It
+follows the range \[`ut.u1_lo`,`ut.u1_up`\] for $\log\varrho$ and
+\[`ut.u2_lo`,`ut.u2_up`\] for $\log\varepsilon$, respectively. Then, the
+2D array `ut.dt2` is allocated and $\{(\log p)_{ij}\}$ data is assigned.
+`ut` is returned and used by the code to finally create the look-up
+table `_TABLP`. Actually, `_TABLP` is a (global) pointer of struct type
+`TAB` (declared in the header file `tabular.h`) and represents the
+tabulated EOS $\log p(\log\varrho,\log\varepsilon)$. Correspondingly,
+coding `eosTabTemperatureUser()` for $\log T$ would provide the input
+for the look-up table `_TABLT`. `_TABLP` and `_TABLT` are used in the
+code function `TDeos()` in module `utilTD.c`.
+
+A tabulated EOS is enabled by appropriate choice of the parameter
+`_C.eos` in the file `nirvana.par` below category PHYSICS
+SPECIFICATIONS.
+
+### User-defined initial mesh refinement and refinement control
+
+For certain problems it may be an advantage to start a simulation with a
+pre-refined mesh in parts of the computational domain, or to have some
+control over mesh refinement in certain domain parts during runtime.
+This can be done with help of the functions `initDomainUser()` and
+`checkDomainUser()` in modules `initDomainUser.c` and
 `checkDomainUser.c`.
 
 #### Initial mesh refinement
 
-The module `initDomainUser.c` serves as template to initialize a
-user-defined refined mesh at simulation start. `initDomainUser()` is
-called by the mesh refinement algorithm. Function arguments are the
-coordinates of a testpoint (`x`,`y`,`z`):
-
-      level=initDomainUser(x,y,z);
-
-The user must return the mesh refinement level (`level`) with which the
-vicinity of the given testpoint should be resolved with the limitation
-\|`level`\|≤`_C.imr` where code paramter `_C.imr` is
-the maximum allowed initial mesh refinement level specified by the user
-in the parameter interface `nirvana.par` under the category
-MESH REFINEMENT. In practice, the user defines spatial domains in
-`initDomainUser.c` via mathematical relations, checks whether the
-testpoint is contained and then returns the requested refinement level.
-The following code fragment illustrates the procedure for a spherical
-interface at radius 0.1 relative to the center
-(*x*, *y*, *z*)=(0.3, 0, 0) to be resolved with refinement level 4:
+In order to refine the mesh initially in a certain spatial domain the
+user must code the function
+
+      level=initDomainUser(x,y,z)
+
+taking arguments (`x`,`y`,`z`) which are coordinates of a testpoint. The
+user must check whether the testpoint lies within its specified domain
+and return the requested refinement level `level` for that domain.
+
+`level` is limited by the constraint
+$|\mathtt{level}|\le\mathtt{\_C.imr}$ where `_C.imr` is the maximum
+allowed initial mesh refinement level as specified in file `nirvana.par`
+below category MESH REFINEMENT. In practice, the user defines spatial
+domains in `initDomainUser()` via mathematical relations. Here is an
+example which illustrates the procedure. It refines the vicinity of a
+spherical interface with radius $0.1$ centered at $(x,y,z)=(0.3,0,0)$
+with refinement level 4:
 
     int initDomainUser(double x, double y, double z)
     /**
-     * Template for the specification of user-defined spatial domains for the 
-     * purpose of initial mesh refinement.
+     * Template for the specification of user-defined spatial domains for initial
+     * mesh refinement.
      *
      * @param x,y,z  spatial coords of testpoint
      *
@@ -1773,14 +1792,15 @@ interface at radius 0.1 relative to the center
      *
      * @author  Udo Ziegler, AIP
      *
-     * @note  Domains must be defined via mathematical relations. The function
-     *        then checks whether a testpoint (x,y,z) lies inside a specified 
-     *        domain. If this is the case the user must return the requested 
-     *        refinement level for this domain. Negative values (level<0) 
-     *        indicate that the domain should be statically refined (i.e. there
-     *        is no derefinement during runtime). For positive values the refined
-     *        domain is subject to subsequent derefinement according to the regular 
-     *        criteria. Multiple domains are allowed.
+     * @note  Domains within the computational grid are to be defined by the user 
+     *        via mathematical relations. This function simply checks then whether 
+     *        a testpoint (x,y,z) lies inside the specified domain(s). If yes, the 
+     *        user must return the requested refinement level for the domain. 
+     *        Especially, if negative values (level<0) are returned the domain is
+     *        assumed to be permanently refined with |level| (i.e., there
+     *        is no derefinement during runtime). Otherwise, if level>0, the 
+     *        initially refined domain is subject to derefinement during runtime.
+     *        Multiple domains are possible.
      */
     {
       int level=0;
@@ -1798,96 +1818,107 @@ interface at radius 0.1 relative to the center
     return level;
     }
 
-The coding here realizes a buffer zone approach by placing three hollow
-globes of decreasing thickness around the interface. This ensures that
-the interface is safely resolved with the requested degree of
-refinement.
-
-*Note: The use of buffer zones is strongly recommended in the setup of
-highly refined, localized regions since the refinement algorithm does
-not allow jumps in refinement more than one level.*
-
-Further examples can be found in the testproblem
-`/nirvana/testproblems/MHD/problem16` and, for a more complicated
-situation, in the testproblem `/nirvana/testproblems/GRAVITY/problem1`.
-
-User-specified refined domains are normally subject to subsequent
-derefinement according to the standard refinement criteria. However, if
-the returned refinement level is given a negative sign
-(`level`&lt;0) it tells the calling function that the
-user-specified domain should be regarded as statically refined. This
-implies that there is no subsequent derefinement of this domain during
-runtime.
-
-#### Restricted mesh refinement
-
-The module `checkDomainUser.c` serves as template for a user-defined
-control mechanism which allows to restrict mesh refinement during
-runtime. `checkDomainUser()` is called by the mesh refinement algorithm.
-Function arguments are the coordinates of a testpoint (`x`,`y`,`z`):
-
-      level=checkDomainUser(x,y,z);
-
-The user must return the refinement control parameter `level`
-(initialized to 0) for the vicinity of the testpoint (`x`,`y`,`z`). A
-return value `level`=0 has no effect on mesh refinement. A
-return value `level`=−1 prohibits mesh refinement. A return
-value `level`&gt;0 forces mesh refinement up to refinement
-level `level` with the limitation `level`≤`_C.amr`
-where code parameter `_C.amr` is the maximum allowed mesh refinement
-level specified by the user in the parameter interface `nirvana.par`
-under category MESH REFINEMENT. Like in the `initDomainUser.c`
-interface the user can define spatial domains via mathematical
-relations, checks whether the testpoint is contained and then returns
-the requested refinement control parameter.
+Note that a buffer zone is realized by placing three spherical shells of
+decreasing thickness around the interface. This is needed because it
+ensures that the interface can safely be resolved with the requested
+degree of refinement level. Further examples can be found in the
+testproblem `/nirvana/testproblems/MHD/problem16` and, for a more
+complicated situation, in the testproblem
+`/nirvana/testproblems/GRAVITY/problem1`.
 
-An example can be found in the testproblem
-`/nirvana/testproblems/GRAVITY/problem3`.
+**Important:** Working with buffer zones is strongly recommended in the
+setup of highly refined, localized regions since the refinement
+algorithm does not allow jumps in refinement across two level. A
+sufficient buffer zone thus ensures that the refinement hierarchy can be
+established smoothly.
+
+An initially refined domain is, by default, subject to derefinement
+during runtime according to the standard refinement criteria. However,
+it can also be marked as static by returning a negative value
+$\mathtt{level}<0$ with $|\mathtt{level}|$ then the requested refinement
+level. No derefinement occurs during runtime but, on the other hand,
+further dynamic refinement on the initially refined domains is possible.
+
+#### Mesh refinement control
 
-## User-defined background magnetic field
+The user can control mesh refinement in spatial domains during runtime
+by the function
 
-The module `sourceB0User.c` serves as template
-for coding a background magnetic field, **B**<sub>0</sub>, in the context
-of the magnetic field splitting formalism.
-**B**<sub>0</sub> must be a time-independent, divergence-free potential field, i.e.,
+      level=checkDomainUser(x,y,z)
 
-d**B**<sub>0</sub>/dt=0,
+taking arguments (`x`,`y`,`z`) which are the coordinates of a testpoint
+as in `initDomainUser()`. As in `initDomainUser()` the user defines
+spatial domains via mathematical relations, checks whether the testpoint
+lies inside a certain domain or not, and returns the corresponding
+control value `level` for that domain:
 
-∇⋅**B**<sub>0</sub>=0
+-   $\mathtt{level}=0$: no restriction for mesh refinement
 
-∇×**B**<sub>0</sub>=0.
+-   $\mathtt{level}=-1$: prohibits mesh refinement on that domain
+
+-   $\mathtt{level}>0$: forces step-by-step mesh refinement on that
+    domain up to level `level` subject to the constraint
+    $\mathtt{level}\le \mathtt{\_C.amr}$ where parameter `_C.amr` is the
+    maximum allowed mesh refinement level as specified in file
+    `nirvana.par` below category MESH REFINEMENT.
+
+An example can be found in the testproblem
+`/nirvana/testproblems/GRAVITY/problem3`.
 
+### User-defined background magnetic field
+
+The module `sourceB0User.c` containing the function `sourceB0User()`
+serves as template for defining a background magnetic field,
+$\mathbf{B}_0$, in the context of the magnetic field splitting formalism
+(cf. [physics
+guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
+$\mathbf{B}_0$ must be a time-independent, divergence-free potential
+field, i.e., 
+$$\partial_t\mathbf{B}_0=0$$,
+$$\nabla\!\cdot\!{\mathbf{B}_0}=0$$ and
+$$\nabla\!\times\!{\mathbf{B}_0}=0$$. 
 In the call of `sourceB0User()`
-the function arguments are the pointer to the array `B0` of magnetic field
-components and the coordinates `x,y,z`:
+the function arguments are the pointer to the array `B0` of magnetic
+field The function
 
-      sourceB0User(B0,x,y,z);
+      sourceB0User(B0,x,y,z)
 
-The user must define the components
-`B0[0]`, `B0[1]` and `B0[2]` in the x-,y- and z-direction, respectively,
-as a function of (x,y,z).
+takes arguments `B0`, a pointer to a 3-element vector for the components
+of the magnetic field, and the spatial coordinates `x,y,z` of a point.
+The user must define the components `B0[0]`=$B0_x$, `B0[1]`=$B0_y$ and
+`B0[2]`=$B0_z$, respectively,
 
-An example definition for a background field (magnetic dipole) can be found in
-testproblem `/nirvana/testproblems/MHD/problem30`.
+An example definition for a background field (a magnetic dipole field)
+can be found in testproblem `/nirvana/testproblems/MHD/problem30`.
 
-The B-splitting scheme is enabled
-by setting the macro `B_FIELD_SPLITTING=YES` in the user interface
-[nirvanaUser.h](3.2-User-interfaces#user-controllable-macros).
+The B-splitting scheme is enabled by setting the macro
+`B_FIELD_SPLITTING=YES` in the header file `nirvanaUser.h` (cf.
+[User-controllable macros](#user-controllable-macros)).
 
-## Specification of NCCM parameters
+### User-defined MHD driver
 
-The parameter file `NCCM.par` serves as user interface to the
-multi-species framework/NCCM. `NCCM.par` is grouped into the category
-sections
+The function `driverMHDUser()` in module `driverMHDUser.c` is thought to
+allow a modification of the MHD flow during the integration cycle.
+Typically, this is a flow driver, or to realize special conditions
+inside the active domain.
 
--   **SPECIES SPECIFICATION** – selection of species.
+An example is given in testproblem
+`/nirvana/testproblems/MHD/problem26`. Here, `driverMHDUser()`
+implements a time-dependent velocity driver to initialize the
+propagation of Alfven wings.
 
--   **REACTIONS SPECIFICATION** – selection of chemical reactions.
+### Specification of NCCM parameters
 
--   **THERMAL PROCESSES SPECIFICATION** – selection of microphysical
-    thermal processes.
+The parameter file `NCCM.par` serves as user interface to NIRVANA's
+multi-species framework and chemistry and cooling model (NCCM).
+`NCCM.par` is grouped into the sections
 
-A section starts with the category name preceded by the ’`>`’-character.
+-   **SPECIES SPECIFICATION** -- selection of species.
+
+-   **REACTIONS SPECIFICATION** -- selection of chemical reactions.
+
+-   **THERMAL PROCESSES SPECIFICATION** -- selection of microphysical
+    thermal processes.
 
 **SPECIES SPECIFICATION**:
 
@@ -1902,40 +1933,39 @@ A section starts with the category name preceded by the ’`>`’-character.
             .          .
             .          .
 
-The following species are currently contained in the multi-species
-framework:
+In the SPECIES SPECIFICATION section the user can select individual
+species. The following species are currently contained:
 
 -   chemical elements H, D, He, C, N, O, Ne, Mg, Si, Fe and their ions
 
--   negatively charged elements H<sup>−</sup>, C<sup>−</sup>,
-    O<sup>−</sup>
-
--   molecules H<sub>2</sub>, H<sub>2</sub>+, H<sub>3</sub><sup>+</sup>,
-    HD, O<sub>2</sub>, C<sub>2</sub>, O<sub>2</sub><sup>+</sup>, CH,
-    CH<sup>+</sup>, CH<sub>2</sub>, CH<sub>2</sub><sup>+</sup>,
-    CH<sub>3</sub><sup>+</sup>, OH, OH<sup>+</sup>, CO, CO<sup>+</sup>,
-    H<sub>2</sub>O, H<sub>2</sub>O<sup>+</sup>,
-    H<sub>3</sub>O<sup>+</sup>, HCO<sup>+</sup> and HOC<sup>+</sup>
-
-Each specie is defined by a unique label (`LABEL`), its chemical
-signature (H\_\_D\_\_He\_C\_\_N\_\_O\_\_Ne\_Mg\_Si\_Fe), i.e., the
-chemical elements a specie consists of, and its charge (`CHARGE`) in
-units of the elementary charge e. These basic species properties should
-*NOT* be changed here by the user. A specie is selected (unselected) by
-setting a `+` (`-`) sign in front of the specie label. Optionally, the
-user can assign each specie its initial mass fraction (`MASSFRAC`) which
-considerably simplifies the definition of IC for the species number
-densities in the user interface `configUser.c` (cf.
-\[Defining IC\](3.2-User-interfaces\#defining-initial-conditions)). The
-species properties are collected in a substructure `_C.X[is]` of the
-public structure `_C` where `is` denotes the species index.
-
-*Note 1: If the sum of `MASSFRAC` over all species does not amount to
-unity the numbers are proportionally rescaled. This changes their
-values.*
-
-*Note 2: A specie electron with label e- is automatically generated by
-the code.*
+-   negatively charged elements H$^-$, C$^-$, O$^-$
+
+-   molecules H$_2$, H$_2$+, H$_3^+$, HD, O$_2$, C$_2$, O$_2^+$, CH,
+    CH$^+$, CH$_2$, CH$_2^+$, CH$_3^+$, OH, OH$^+$, CO, CO$^+$, H$_2$O,
+    H$_2$O$^+$, H$_3$O$^+$, HCO$^+$ and HOC$^+$
+
+Each species is defined by a unique label `LABEL` (eg., He+, Fe+16,
+HOC+), its chemical signature (H\_\_D\_\_He_C\_\_N\_\_O\_\_Ne_Mg_Si_Fe),
+i.e., the chemical elements is made of, its charge `CHARGE` in units of
+the elementary charge e, and its initial mass fraction `MASSFRAC`.
+`LABEL`, the chemical signature and CHARGE should **NOT** be changed by
+the user. Species are selected (unselected) by setting a `+` (`-`) sign
+in front of the species label. Optionally, the user can assign each
+species its initial mass fraction which considerably simplifies the
+definition of IC for the species number densities in `configUser.c` (cf.
+[Defining initial conditions](#defining-initial-conditions)). The sum of
+mass fractions should be 1. If not, the mass fractions are automatically
+rescaled changing the given value. The species properties are collected
+in a substructure `_C.X[is]` below the master control structure `_C`
+where `is` denotes the species index. Within the specified procedure, it
+is possible to add further species.
+
+**Important:** A species 'electron' with label e- is automatically
+generated by the code and not shown in `NCCM.par`.
+
+If the multi-species framework is enabled (section PHYSICS
+SPECIFICATIONS in file `nirvana.par`) an advection(-reaction) equation
+is solved for each selected species.
 
 **REACTIONS SPECIFICATION**:
 
@@ -1951,12 +1981,13 @@ the code.*
             .          .
             .          .
 
-The chemical reactions in section REACTIONS SPECIFICATION are grouped
-into the classes:
+In the REACTIONS SPECIFICATIONS section chemical reactions between
+species are selected. The predefined chemical reactions so far are
+grouped into the following categories:
 
 -   subnetworks for the ionization structure of each element
 
--   subnetworks for H<sub>2</sub> formation and HD formation
+-   subnetworks for H$_2$ formation and HD formation
 
 -   C-bearing and O-bearing chemical cycles
 
@@ -1964,19 +1995,24 @@ into the classes:
 
 -   dust-assisted reactions
 
+Please consult the [chemistry
+guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/ChemistryGuide.pdf)
+for more details on the current chemical network of NIRVANA.
+
 Each reaction is characterized by an ID (a running number) followed by
 its definition according to the pattern
 
            STCOEFF REACTANT [& STCOEFF REACTANT] > STCOEFF PRODUCT [& STCOEFF PRODUCT]
 
-with stoichiometric coefficient STCOEFF, reactants REACTANT and products
-PRODUCT. A reaction is selected (unselected) by setting a `+ (-)` sign
-in front of the reaction ID.
-
-*Note: The code parser checks each reaction for charge conservation and
-mass conservation as well as the existence of the labels of involved
-species. This rection check can be disabled by setting the macro
-`CHECK_REACTIONS` to `NO` in `nirvanaUser.h`.*
+with stoichiometric coefficients STCOEFF, reactants REACTANT and
+products PRODUCT. A reaction is selected (unselected) by setting a
+`+ (-)` sign in front of the reaction ID. Within the specified
+nomenclature it is possible to add further reactions. Note that the code
+parser checks each reaction for charge conservation and mass
+conservation, as well as the existence of the involved species (by
+checking for `LABEL` in the SPECIES SPECIFICATIONS section). Reaction
+checking can be disabled by setting the macro `CHECK_REACTIONS=NO` in
+header file `nirvanaUser.h`.
 
 **THERMAL PROCESSES SPECIFICATION**:
 
@@ -2000,17 +2036,18 @@ species. This rection check can be disabled by setting the macro
      N   >H2O: CE COOLING
      N   >OH: CE COOLING
 
-The NCCM currently contains the following microphysical thermal
-processes:
+In the THERMAL PROCESSES SPECIFICATION section the user can select
+individual thermal process associated with the chemical network. The
+NCCM currently includes the following processes:
 
--   Bremsstrahlung for the H, He, O, Ne and Fe ions
+-   Bremsstrahlung for H, He and He ions
 
 -   chemical- and atomic line cooling of H, D and He
 
 -   fine-structure line cooling of the metals C, N, O, Ne, Mg, Si and Fe
 
 -   chemical-, rotovibrational line- and collision-induced cooling of
-    H<sub>2</sub>
+    H$_2$
 
 -   rotovibrational line cooling of HD
 
@@ -2018,36 +2055,49 @@ processes:
 
 -   rotovibrational line cooling of CO
 
--   rotovibrational line cooling of H<sub>2</sub>O
+-   rotovibrational line cooling of H$_2$O
 
 -   rotational line cooling of OH
 
-A thermal process is activated (deactivated) by specifying the value Y
-(N). Every text in a line following the ’`>`’-character is ignored by
-the parser and serves for comments.
+A thermal process is selected (unselected) by specifying the value Y
+(N).
 
-## User-controllable macros
+### User-controllable macros
 
-Besides the parameters collected in the user interfaces `nirvana.par`
-and `NCCM.par` there are a number of further parameters which can be
-controlled by the user. In practice, however, the user rarely has to
-change them. These parameters are defined in the header file
-`nirvanaUser.h` as macros and it are grouped into the categories:
+The files `nirvana.par` and `NCCM.par` serve for the specification of
+primary simulation parameters. controlled by the user. In practice,
+however, the user rarely has to change them. There are further
+controllable parameters defined as macros in the header file
+`nirvanaUser.h`. It are grouped into the categories/subcategories:
 
 -   NUMERICS-RELATED MACROS
 
+    -   FV SCHEME
+
+    -   MESH REFINEMENT
+
+    -   GRAVITY/MULTIGRID
+
+    -   RKL SCHEME
+
+    -   HEATLOSS SOLVER
+
+    -   CHEMICAL REACTIONS SOLVER
+
 -   PHYSICS-RELATED MACROS
 
 -   NCCM-RELATED MACROS
 
--   MACROS RELATED TO ANALYTICAL EOS
+-   ANALYTICAL EOS-RELATED MACROS
 
 The following table compiles a complete list of macros including a short
-description of their meaning and, in some cases, their permissable
+description of their meaning and, in some cases, their allowed range of
 values (given in brackets).
 
 **NUMERICS-RELATED MACROS**:
 
+-   FV SCHEME
+
     -   `SPACE_ORDER` (2): spatial order of the numerical scheme
 
     -   `LIM(a)` ({MM(a),MC(a),VL(a)}): type of slope limiter in the
@@ -2058,65 +2108,111 @@ values (given in brackets).
 
     -   `PSI` (\[0.5,1\]): parameter in multi-d limiter
 
--   `MAXLEVEL` (&lt;128): maximum refinement level
+    -   `HLLD_PRESSURE_CORR` ({YES,NO}): option to enable/disable the
+        low Mach pressure correction term in the HLLD Riemann solver
+        according to [@MM21]
+
+    -   `BORIS_CORR` ({YES,NO,AUTO}): option to enable/disable the Boris
+        correction (cf. [Code features](#code-features))
+
+    -   `BORIS_CRED`: reduced speed of light in Boris correction
+
+    -   `BORIS_CRED_MIN`: lower limit of reduced speed of light in units
+        $2|{\bf v}|+c_s$
+
+    -   `BORIS_AUTO_CA_MAX`: Alfven speed threshold in units
+        $2|{\bf v}|+c_s$ above which the Boris correction is avtivated
+        when `BORIS_CORR`=`AUTO`
 
--   `MG_ITMAX`: maximum number of iterations allowed in the multigrid
+-   MESH REFINEMENT
+
+    -   `MAXLEVEL` ($<128$): maximum refinement level
+
+    -   `NOBR` ({YES,NO}): option to prevent mesh refinement at domain
+        boundaries
+
+-   GRAVITY/MULTIGRID
+
+    -   `GR_RES_NORM` (0,1,2): norm to be used in calculating the
+        residual in the multigrid solver:
+
+        -   0: 1-norm (normalized to initial 1-norm)
+
+        -   1: relative 1-norm
+
+        -   2: relative max-norm
+
+    -   `MG_ITMAX`: maximum number of multigrid iterations before
+        termination
+
+    -   `MG_TOL` (typical: $10^{-6}$): error tolerance in the multigrid
         solver
 
--   `MG_TYPE` (MULT): used type of multigrid solver
+    -   `MG_N_SUBLEVELS` (1,2): number of sublevels (below zero level)
+        to be used in the multigrid solver
 
--   `MG_TOL` (typical: 10<sup>−6</sup>): error tolerance in the
-    multigrid solver
+    -   `MG_SMOOTHER`
+        ({AUTO,POINT,XLINE,YLINE,ZLINE,XYPLANE,XZPLANE,YZPLANE}):
+        smoother type
 
--   `RKL_COURANT_EXPL` (&lt;0.5): explicit Courant number used in the
-    RKL solver
+    -   `MG_ORDER` ({LEX,RB}): relaxation order
 
--   `RKL_MAX_COURANT` (typical: &lt;1000): maximum Courant-like number
-    allowed in the RKL solver
+    -   `MG_BASE_ITMAX`: base level solver maximum iterations
 
--   `RKL_DT_LIM`: RKL-related timestep limit as fraction of the
-    dynamical timestep
+    -   `MG_BASE_RES`: base level residual threshold
 
--   `HEATLOSS_TOL` (typical: 10<sup>−5</sup>): relative error
-    tolerance in heatloss solver
+    -   `MG_BASE_ORDER` ({LEX,RB}): base level relaxation order
 
--   `HEATLOSS_ATOL`: absolute error tolerance in heatloss solver
+    -   `MG_BASE_SOR` ({YES,NO}): option to enable overrelaxation in the
+        base level solver
 
--   `REACTIONS_TOL` (typical: 10<sup>−5</sup>): relative error
-    tolerance in NCCM solver
+    -   `MG_BASE_CHEB_ACC` ({YES,NO}): option to enable Chebyshev
+        acceleration in the
 
--   `REACTIONS_ATOL_X`: absolute error tolerance for number densities in
-    NCCM solver
+    -   `MG_MAX_ANISO_X` ($\ge 1$): allowed y,z-to-x anisotropy
 
--   `REACTIONS_ATOL_T`: absolute error tolerance for the temperature in
-    NCCM solver
+    -   `MG_MAX_ANISO_Y` ($\ge 1$): allowed z-to-y anisotropy base level
+        solver
 
--   `HLLD_PRESSURE_CORR` ({YES,NO}): option to enable/disable the low Mach 
-    pressure correction term in the HLLD Riemann solver according to 
-    \[[MM21](3.2-User-interfaces#references)\]
+-   RKL SCHEME
 
--   `BORIS_CORR` ({YES,NO,AUTO}): option to enable/disable the Boris correction
-    (cf. [Functionality overview](3.1-Code-basics#functionality-overview)) 
+    -   `RKL_COURANT_EXPL` ($<0.5$): explicit Courant number used in the
+        RKL solver
 
--   `BORIS_CRED`: reduced speed of light in Boris correction
+    -   `RKL_MAX_COURANT` (typical: $<1000$): maximum Courant-like
+        number allowed in the RKL solver
 
--   `BORIS_CRED_MIN`: lower limit of reduced speed of light in units 
-    2|**v**|+c<sub>s</sub> in Boris correction
+    -   `RKL_DT_LIM`: RKL-related timestep limit as fraction of the
+        dynamical timestep
+
+-   HEATLOSS SOLVER
+
+    -   `HEATLOSS_TOL` (typical: $10^{-5}$): relative error tolerance in
+        heatloss solver
+
+    -   `HEATLOSS_ATOL`: absolute error tolerance in heatloss solver
+
+-   CHEMICAL REACTIONS SOLVER
+
+    -   `REACTIONS_TOL` (typical: $10^{-5}$): relative error tolerance
+        in NCCM solver
+
+    -   `REACTIONS_ATOL_X`: absolute error tolerance for number
+        densities in NCCM solver
 
--   `BORIS_AUTO_CA_MAX`: Alfven speed threshold in units 2|**v**|+c<sub>s</sub> 
-    above which the Boris correction is avtivated when `BORIS_CORR=AUTO`
+    -   `REACTIONS_ATOL_T`: absolute error tolerance for the temperature
+        in NCCM solver
 
 **PHYSICS-RELATED MACROS**
 
--   `CENTRIFUGAL_FORCE` ({YES,NO}): option to enable/disable the
-    centrifugal force term
+-   `CENTRIFUGAL_FORCE` ({YES,NO}): option to include the centrifugal
+    force term
 
 -   `COND_FORCE_ISO` ({YES,NO}): option to force isotropy in thermal
     conduction
 
--   `B_FIELD_SPLITTING` ({YES,NO}): option to enable/disable the magnetic 
-    field splitting formalism
-    (cf. [Functionality overview](3.1-Code-basics#functionality-overview)) 
+-   `B_FIELD_SPLITTING` ({YES,NO}): option to use the magnetic field
+    splitting formalism (cf. [Code features](#code-features))
 
 **NCCM-RELATED MACROS**
 
@@ -2126,28 +2222,28 @@ values (given in brackets).
 
 -   `ADUST`: dust abundance relative to Milky Way
 
--   `ISRF_G` (Milky Way: 1.13): Habing strength of interstellar
-    radiation field
+-   `ISRF_G`: Habing strength of interstellar radiation field (1.13)
 
 -   `DBNCCM` ({YES,NO}): NCCM database availability flag
 
 -   `CHECK_REACTIONS` ({YES,NO}): option to check chemical reactions
     with regard to mass conservation, charge conservation and existence
-    of species labels
+    of involved species names
 
-**MACROS RELATED TO ANALYTICAL EOS**
+**ANALYTICAL EOS-RELATED MACROS**
 
--   `PUSR(rho,eth)`: analytic expression for the pressure *p* as a
-    function of 𝜚 (`rho`) and *ε* (`eth`)
+-   `PUSR(rho,eth)`: user-defined expression for the pressure $p$ as a
+    function of $\varrho$ (`rho`) and $\varepsilon$ (`eth`)
 
--   `CS2USR(rho,eth)`: analytic expression for the square of sound speed
-    *c*<sub>*s*</sub><sup>2</sup>(𝜚, *ε*)
+-   `CS2USR(rho,eth)`: user-defined expression for the square of sound
+    speed $c_s^2(\varrho,\varepsilon)$
 
--   `TUSR(rho,eth)`: analytic expression for the temperature *T*(𝜚, *ε*)
+-   `TUSR(rho,eth)`: user-defined expression for the temperature
+    $T(\varrho,\varepsilon)$
 
--   `ETUSR(rho,eth)`: analytic expression of the thermal energy density
-    *ε*(𝜚) as a function of 𝜚 in cases without energy equation.
-    Otherwise identity.  
+-   `ETUSR(rho,eth)`: user-defined expression of the thermal energy
+    density $\varepsilon(\varrho)$ as a function of $\varrho$ in
+    simulations without energy equation. Otherwise identity.
 
 ### References