epoch_dev.tex

\include{epoch_head}
\externaldocument{epoch_user}

\begin{document}
\includepdf{images/title_page_dev}
{
  \fontfamily{phv}\selectfont\input{epoch_dev_title}
}
\fontfamily{garamond}\selectfont%
\tableofcontents%
\newpage%
\DefineShortVerb{\#}
\fvset{formatcom=\color{warwickred}}

\section{{\EPOCH} developers manual}

This section of the manual is aimed at people who intend to edit the {\EPOCH}
source code to extend or modify existing features, add new diagnostics or
develop new physics packages. It is expected that anyone reading this part of
the manual will be familiar with the material covered in the main manual,
particularly the notices about the layout of particle, particle\_list and
particle\_species structures in \sect{partrep} of the user manual.

A quick note: Some files have the normal Fortran file extension .f90, while
some have the slightly unusual .F90. The difference is that files with the
.F90 extension are passed through the preprocessor before they are compiled
allowing the use of precompiler directives (the \#ifdef commands).

\scaledcapimage{./images/coreblock}{coreblock}{Block diagram of the core
  routines in {\EPOCH}}{0.5}

\section{General layout of the {\EPOCH} code}

The names of the source files in {\EPOCH} are fairly self explanatory but, for
clarity, they are explained here.

\subsection{Directories}
All source files are contained in the \inlinecode{src} directory and its
subdirectories. There is a stylistic reason for the layout of the files, which
is explained here

\begin{itemize}
\item \inlinecode{src} - Files in this directory are the core files for the
  basic {\EPOCH} code, such as the field solvers, the particle pusher, the
  boundary conditions and the lasers.
\item \inlinecode{src/deck} - Files in this directory are responsible for
  dealing with the permanent input deck parser and include the core parts of
  the deck handler, and also the routines which deal with the blocks in the
  input deck files.
\item \inlinecode{src/housekeeping} - Files in this directory deal with those
  parts of the code operation which are not physics; including the load
  balancer, the MPI setup routines and the moving window.
\item \inlinecode{src/include} - This directory contains the shape function
  code fragments which are inserted into the particle push at compile time.
\item \inlinecode{src/io} - The files involved in all I/O activities, including
  the distribution functions and the particle probes.
\item \inlinecode{src/parser} - The files for the maths expression parser are
  in this directory, including both the core implementation of the shunting yard
  algorithm and the routines for implementing the permanent functions,
  constants and operators for the input deck.
\item \inlinecode{src/physics\_packages} - Contains routines which implement
  additional physics for the code.
\item \inlinecode{src/user\_interaction} - Contains any Fortran routines which
  a user has to modify to use the code with internal initial conditions, or to
  temporarily extend the maths parser or the input deck.
\end{itemize}


\subsection{The files in \inlinecode{src}}
\begin{itemize}
\item\inlinecode{boundary.f90} - Includes all boundary conditions except laser
  and transmissive boundaries; including field and particle MPI boundaries, and
  field and particle domain boundaries.
\item\inlinecode{epoch\{n\}d.F90} - Main driver for the code. Reading this
  routine gives the basic layout of the code flow.
\item\inlinecode{fields.f90} - The Maxwell field solver.
\item\inlinecode{laser.f90} - Includes laser and transmissive boundary
  conditions for each boundary and also the housekeeping routines for the laser
  objects.
\item\inlinecode{particles.F90} - The particle pusher.
\item\inlinecode{shared\_data.F90} - This file includes all the global variable
  and type definitions. Usually new variables should be defined in this file.
\item\inlinecode{gen\_commit\_string} - This is a script used to generate an ID
  string when compiling the code.
\item\inlinecode{gen\_src\_module} - This is a script which is used at build
  time to generate a Fortran module containing the source code. This is used for
  embedding the {\EPOCH} source code into restart dumps.
\end{itemize}

\subsection{The files in \inlinecode{src/deck}}
\label{sec:src_deck}
\begin{itemize}
\item\inlinecode{deck.F90} - The main input deck routines. Deals with opening
  files, reading data and MPI distribution of the data to all processes. Also
  includes the routines which deal with calling the right reader routines to
  deal with a given block.
\item\inlinecode{deck\_boundaries\_block.f90} - Reader routine for the
  ``boundaries'' block of the input deck.
\item\inlinecode{deck\_constant\_block.f90} - Reader routine for ``constants''
  blocks in the input deck.
\item\inlinecode{deck\_control\_block.f90} - Reader routine for ``control''
  block in the input deck.
\item\inlinecode{deck\_dist\_fn\_block.f90} - Reader routine for ``dist\_fn''
  blocks in the input deck.
\item\inlinecode{deck\_fields\_block.f90} - Reader routine for ``fields''
  blocks in the input deck.
\item\inlinecode{deck\_io\_block.F90} - Reader routine for the ``io'' block in
  the input deck.
\item\inlinecode{deck\_laser\_block.f90} - Reader routine for ``laser'' blocks
  in the input deck.
\item\inlinecode{deck\_particle\_probe\_block.F90} - Reader routine for
  ``particle\_probe'' blocks in the input deck.
\item\inlinecode{deck\_species\_block.F90} - Reader routine for ``species''
  blocks in the input deck.
\item\inlinecode{deck\_window\_block.f90} - Reader routine for the ``window''
  block in the input deck.
\item\inlinecode{strings.f90} - Basic string handling routines such as
  ``str\_cmp'' and routines for converting strings to numbers WITHOUT using the
  maths parser are covered in this routine.
\item\inlinecode{strings\_advanced.f90} - The routines which pass maths along
  to the maths parser routines are here.
\end{itemize}

\subsection{The files in \inlinecode{src/housekeeping}}
\begin{itemize}
\item\inlinecode{balance.F90} - Contains the routines for the load balancer
  and related routines.
\item\inlinecode{current\_smooth.F90} - Contains the current smoothing routines.
\item\inlinecode{dummy\_encoded\_source.f90} - A dummy module used when the
  code is being built without the ability to write the source code into restart
  dumps.
\item\inlinecode{mpi\_routines.F90} - Contains the routines dealing with the
  setup of the MPI layer and the creation of the communicator. Also allocates
  all arrays for the first time before load balancing.
\item\inlinecode{mpi\_subtype\_control.f90} - Contains the routines that setup
  the mpi types required by the I/O subsystem.
\item\inlinecode{particle\_pointer\_advance.f90} - Contains subroutines which
  walk through the lists of particles and species for I/O purposes.
\item\inlinecode{partlist.F90} - Contains the routines which deal with the
  particle lists which are used for inter-processor communication of particles.
\item\inlinecode{random\_generator.f90} - Contains the random number generator
  routines.
\item\inlinecode{setup.F90} - Deals with the setup of the grids and domains and
  restarting from previous output dumps.
\item\inlinecode{shape\_functions.F90} - Contain the particle shape functions
  used for calculating the particle weighting.
\item\inlinecode{split\_particle.F90} - Is the implementation of a
  demonstration of particle splitting routines.
\item\inlinecode{utilities.f90} - Contains growable arrays used by the species
  block parser.
\item\inlinecode{version\_data.F90} - Contains version information for the
  current {\EPOCH} code.
\item\inlinecode{welcome.F90} - The routine which prints the banner message and
  compiler options info.
\item\inlinecode{window.F90} - The routines which deal with the moving window.
\end{itemize}

\subsection{The files in \inlinecode{src/io}}
\begin{itemize}
\item\inlinecode{calc\_df.F90} - Despite the slightly confusing name, this
  subroutine deals with derived functions like number density, charge density
  and mass density.
\item\inlinecode{diagnostics.F90} - Contains the routines which actually dump
  the data, decide what to dump and also the routine to calculate the timestep.
\item\inlinecode{dist\_fn.F90} - Contains the routines to calculate the
  distribution functions, and also the routines handling the requests for
  distribution functions.
\item\inlinecode{iterators.F90} - Contains the iterator functions used to write
  particle data into SDF files.
\item\inlinecode{probes.F90} - Contains the routines which write the data from
  the particle probes. Also includes the routines which deal with user requests
  to add new particle probes.
\item\inlinecode{sdf\_common.f90} - Core part of the SDF file format. Contains
  constants and definitions needed by the SDF format. If attempting to
  implement a SDF reader then look here for values of named constants.
\item\inlinecode{sdf\_control.f90} - Core part of the SDF file format. Contains
  the routines needed to open and close files etc.
\item\inlinecode{sdf.f90} - Module interface for the SDF file format. This file
  specifies the application interface to the SDF library.
\item\inlinecode{sdf\_input.f90} - Core part of the SDF file format. Contains
  the general routines needed for reading data from SDF files.
\item\inlinecode{sdf\_input\_cartesian.f90} - Core part of the SDF file format.
  Contains the routines needed to read Cartesian meshes and variables from SDF
  files.
\item\inlinecode{sdf\_input\_functions.f90} - Core part of the SDF file format.
  Contains the routines needed to traverse the block structure of a SDF file
  without reading any block specific metadata.
\item\inlinecode{sdf\_input\_point.f90} - Core part of the SDF file format.
  Contains the routines needed to read particles meshes and variables from SDF
  files.
\item\inlinecode{sdf\_input\_util.f90} - Core part of the SDF file format.
  Contains general routines for reading the block structure of SDF files.
\item\inlinecode{sdf\_job\_info.f90} - Contains routines for generating a
  unique job-id for each simulation run.
\item\inlinecode{sdf\_output.f90} - Core part of the SDF file format. Contains
  the basic routines needed to write data to a SDF file.
\item\inlinecode{sdf\_output\_cartesian.f90} - Core part of the SDF file
  format. Contains the routines needed to write Cartesian meshes and variables.
\item\inlinecode{sdf\_output\_point.f90} - Core part of the SDF file format.
  Contains the routines needed to write particles meshes and variables.
\item\inlinecode{sdf\_output\_util.f90} - Core part of the SDF file format.
  Contains general routines for writing the block structure of SDF files.
\item\inlinecode{simple\_io.f90} - Contains routines for performing the simple
  binary I/O required by species\_external and fields\_external blocks.
\end{itemize}

\subsection{The files in \inlinecode{src/parser}}
\begin{itemize}
\item\inlinecode{evaluate.f90} - Contains the routines which actually evaluate
  a tokenized expression. The core of this is a simple implementation of an
  RPN calculator.
\item\inlinecode{evaluator\_blocks.f90} - Contains the routines which evaluate
  a given token into a numerical values. Actually implements the functions,
  constants and operators in {\EPOCH}'s maths parser.
\item\inlinecode{shunt.F90} - {\EPOCH}'s implementation of the
  ``shunting yard'' algorithm used to simultaneously tokenize the input and
  convert it from infix notation to RPN.
\item\inlinecode{stack.f90} - Deals with routines for pushing onto and popping
  off stacks.
\item\inlinecode{tokenizer\_blocks.f90} - Deal with converting strings found
  in a string being tokenized into tokens. Essentially a large collection of
  ``str\_cmp'' commands testing a string against a known name.
\end{itemize}

\subsection{The files in \inlinecode{src/physics\_packages}}
\begin{itemize}
\item\inlinecode{TABLES} - Subdirectory containing tables
for QED emissions. 
\item\inlinecode{numerics.f90} - Some additional numerics
routines, primarily for the photons (QED) package.
\item\inlinecode{ionise.F90} - Physics package dealing with field ionisation. 
Field ionisation consists of three distinct regimes; multiphoton in which
ionisation is best described as absorption of multiple photons, tunnelling 
in which deformation of the atomic Coulomb potential is the dominant factor,
and barrier suppression ionisation in which the electric field is strong enough
for an electron to escape classically. It is possible to turn off multiphoton or
barrier suppression ionisation through the input deck. 
\item\inlinecode{collisions.F90} - Physics package handling particle collisions
and collisional ionisation. See 
\url{http://iopscience.iop.org/article/10.1088/0741-3335/57/11/113001/pdf}
for details. 
\item\inlinecode{injectors.F90} - Package dealing with particle injection through
a boundary. Currently this handles drifting and non-drifting Maxwellian and
flux-Maxwellian distributions. 
\item\inlinecode{photons.F90} - Package for some QED effects. 
EPOCH can model QED pair production, synchrotron emission and radiation 
reaction as described in Duclous et al 
\url{http://iopscience.iop.org/article/10.1088/0741-3335/53/1/015009} and 
Ridgers et al. 
\url{https://www.sciencedirect.com/science/article/pii/S0021999113008061?via}
\end{itemize}

\subsection{The files in \inlinecode{src/user\_interaction}}
\begin{itemize}
\item\inlinecode{custom\_deck.f90} - This file is where and end user can
  temporarily extend the input deck. Described in the
  ``\nameref{sec:customising}'' section.
\item\inlinecode{custom\_laser.f90} - The file where an end user specifies
  laser time profiles without using the input deck.
\item\inlinecode{custom\_parser.f90} - The file where an end user can
  temporarily add new functions and constants to the input deck. It is
  described in the ``\nameref{sec:customising}'' section of this manual.
\item\inlinecode{helper.F90} - This file contains all the internal workings of
  the autoloader. This is in user\_interaction for historical reasons, since
  early versions of the code required the end user to modify some parts of the
  functions contained in this file. As the autoloader has increased in
  complexity, this has ceased to be the case, so it is likely that soon this
  file will be move to ``housekeeping''.
\item\inlinecode{ic\_module.f90} - This file is where the internal and manual
  initial conditions are set.
\item\inlinecode{particle\_temperature.F90} - Contains the routines for
  thermally loading a particle species.
\end{itemize}

\section{{\EPOCH} makefile}

The makefile supplied with {\EPOCH} is a standard GNU make makefile, which must
be user modified to allow a developer to add new files to the code. {\EPOCH}'s
makefile is quite large, so an explanation of how to add new files and new
directories is given below.

\subsection{Adding a new file to be compiled with {\EPOCH}}
There are three things that must be done to cause {\EPOCH} to compile a new
file and link it into the final code. Assume that you're adding a file called
``newfile.F90''. First, find the line which sets the environment variable
\inlinecode{SRCFILES} and add a new parameter which reads\\
\\
newfile.F90\\
\\ This tells the makefile to compile the final code using your new file, the
next thing to do is to add a line which tells the code about the dependencies
for your file. Lower down in the makefile, you'll find a section with lines
which look like:
\begin{boxverbatim}
balance.o: balance.F90 boundary.o mpi_subtype_control.o partlist.o
\end{boxverbatim}
Add a new line for describing all the FILES (NOT modules) which are used by
your new file. If you USE shared\_data, mpi\_subtype\_control and stack in
your file then the line would look like:
\begin{boxverbatim}
newfile.o: newfile.F90 shared_data.o mpi_subtype_control.o stack.o
\end{boxverbatim}
Note the structure of the line with ONLY the source file for the new file
specified, all other used files specify the intermediate .o files. The
remaining element of the makefile which needs to be modified is to add your
new file as a dependency to all the files which USE modules contained in your
new file. This is achieved very simply by adding ``newfile.o'' to the dependency
list for those files which USE your modules. For example if you've written new
boundary conditions and USE your modules in boundary.f90, you'd just change
the line for boundary.f90 from:
\begin{boxverbatim}
boundary.o: boundary.f90 deck_io_block.o particle_temperature.o partlist.o
\end{boxverbatim}
to
\begin{boxverbatim}
boundary.o: boundary.f90 deck_io_block.o particle_temperature.o partlist.o \
  newfile.o
\end{boxverbatim}

Note that the backslash characters are line continuation marks in makefiles.

\subsection{Adding new directories to {\EPOCH}'s makefile}
If you want to add an entire new directory to the {\EPOCH} compile path then
you need to add it to the definition of the variable \inlinecode{VPATH}.
Remember to use the variable \inlinecode{\$(SRCDIR)} rather than hard-coding
\inlinecode{src} into the path.

\section{{\EPOCH} core programming}

{\EPOCH} is designed so that it can fairly easily be extended while still being
written in (more or less) standard Fortran90 and MPI1.2. This section details
in increasing complexity what a programmer needs to know to extend {\EPOCH} with
new diagnostics, new physics or even new core solvers. The first few entries
in this section range between style points and explanations of fairly trivial
parts of the {\EPOCH} code, but the end of this section gives an overview of how
one would perform major changes to the complete {\EPOCH} core solver.

\subsection{Physical constants}
In order to ensure that different parts of the code run at the same precision
common physical constants are defined in \inlinecode{shared\_data.F90} and any
new physical constants required by extensions to the code should be placed in
the same location. The constants available in the code are
\begin{itemize}
\item pi - Ratio of a circle's circumference to its diameter.
\item q0 - Charge on electron.
\item m0 - Rest mass of electron.
\item c - Speed of light in vacuum.
\item kb - Boltzmann's constant.
\item mu0 - Permeability of free space.
\item epsilon0 - Permittivity of free space.
\item h\_planck - Plank's constant.
\item ev - The value of an electron volt.
\end{itemize}

Any new constants required should be specified in the same place in
\inlinecode{shared\_data.F90}.

\subsection{Important variables, arrays and array length}
As well as the physical constants, there are some important variables which
you will have to use to do any development with {\EPOCH}. As a general note,
since {\EPOCH} is written with separate 1D, 2D and 3D versions, definitions will
be given for the 3D version of the code and irrelevant dimensions should just
be left out.

\subsubsection{Shape and size variables}
\begin{itemize}
\item #INTEGER :: nx, ny, nz# - The number of gridpoints on the current
  processor in each direction. This may change when the load balancer
  activates, so always use these variables rather than local copies.
\item #INTEGER :: nx_global, ny_global, nz_global# - The number of gridpoints
  in each direction of the whole domain. These numbers will never change and
  will be the numbers read in from the input deck.
\item #INTEGER(KIND=8) :: npart_global# - The global number of particles
  specified in the input deck. This is not updated as particles leave the
  domain through boundaries etc. so it is not guaranteed to be accurate.
\item #INTEGER :: n_species# - The number of species of particles specified.
\item #INTEGER :: nsteps# - The maximum number of steps that the core solver
  should take.
\setlength{\emergencystretch}{3em}
\item #INTEGER, DIMENSION(1:nproc{x,y,z}), ALLOCATABLE :: cell_{x,y,z}_min,#
  \linebreak#cell_{x,y,z}_max# - The variables #cell_{x,y,z}_min# and
  #cell_{x,y,z}_max# represent the part of a global array which is held by the
  current processor. Since {\EPOCH} is an MPI code, there doesn't exist a
  single copy of any of the global arrays anywhere, but if there did then each
  processor would be responsible for the slice which runs
  (cell\_x\_min(rank):cell\_x\_max(rank),
  cell\_y\_min(rank):cell\_y\_max(rank))
  in 2D.  These variables are used internally in the load balancer, where it is
  updated, but is also used when calculating distribution functions. Here it is
  used to define the extents of the MPI type which is used to write the
  distribution function to disk.
\end{itemize}

\subsubsection{Input deck variables}
\begin{itemize}
\item #CHARACTER(LEN=entry_length) :: blank# - A special string which the input
  deck parser uses to indicate that it's passing a blank string rather than a
  string read from the deck which just happens to be blank.
\item #INTEGER :: deck_state# - An integer determining the current sweep of
  the input deck by the deck parser.
\item #INTEGER, PARAMETER :: num_vars_to_dump# - A variable describing the
  number of variables which should be selectable in the input deck as possible
  variables to dump.
\item #CHARACTER(LEN=entry_length) :: extended_error_string# - String used by
  some error codes in the deck parser to give more user friendly error
  messages.
\item #INTEGER :: data_dir_max_length# - The maximum number of characters in
  the name of the output directory.
\item #INTEGER :: n_zeros# - The number of leading zeros in the output filenames
  from {\EPOCH}.
\end{itemize}

\subsubsection{Initial conditions (autoloader) variables}
Initial conditions for the autoloader for a given species are described in
{\EPOCH} by the Fortran TYPE \inlinecode{initial\_conditions\_block}. The
definition (in 3D) is:
\begin{boxverbatim}
TYPE initial_condition_block
  REAL(num), DIMENSION(:,:,:), POINTER :: density
  REAL(num), DIMENSION(:,:,:,:), POINTER :: temp
  REAL(num), DIMENSION(:,:,:,:), POINTER :: drift

  REAL(num) :: density_min
  REAL(num) :: density_max
END TYPE initial_condition_block
\end{boxverbatim}

In 2D, the arrays have one fewer index, and in 1D they have two fewer.

\begin{itemize}
\item #REAL(num) :: density# - Number density for the particles in the species.
  When defined runs (-2:nx+3,-2:ny+3,-2:nz+3).
\item #REAL(num) :: temp# - Temperature in Kelvin of the species in space. When
  defined runs (-2:nx+3,-2:ny+3,-2:nz+3,1:3). The final index of the array
  is a direction index, used to give anisotropic thermal distributions.
\item #REAL(num) :: drift# - Velocity drift in $m/s$ of the species in space.
  When defined runs (-2:nx+3,-2:ny+3,-2:nz+3,1:3). The final index of the array
  is the velocity direction component.
\item #density_min# - The minimum density below which the autoloader
  should not load particles.
\item #density_max# - The maximum density above which the autoloader
  should clip the density function.
\end{itemize}

The initial conditions themselves are in the variable
\begin{boxverbatim}
TYPE(initial_condition_block), DIMENSION(:), ALLOCATABLE :: initial_conditions
\end{boxverbatim}
which is allocated to an array of size \inlinecode{1:n\_species}.

\subsubsection{Linked Lists}
Linked lists are a standard computer programming technique which is still
slightly unusual in Fortran, and may well be unfamiliar to many Fortran
programmers. They effectively allow you to have an array of arbitrary length,
although this comes with various trade-offs about memory locality and speed of
accessing elements. The general concept is that of a chain where each link in
the chain only knows about the previous link in the chain and the next link in
the chain. Although there are schemes for doing this in languages which don't
have pointers, the normal method of implementing linked lists is to use
pointers to point to previous and next elements in the list, and this is how
they are implemented in {\EPOCH}. Since both linked lists and Fortran pointers
are slightly esoteric concepts, while being key to the operation of {\EPOCH} a
brief overview of them is presented here.\\

The simplest possible form of a linked list element would be a TYPE which
looks like:
\begin{boxverbatim}
TYPE linked_list
  TYPE(linked_list), POINTER :: next
  TYPE(linked_list), POINTER :: prev
END TYPE linked_list
\end{boxverbatim}

You also have to have a pointer to the start of the list, and to speed up
adding new elements to the list, you normally also keep a pointer to the
last element of the list. Therefore, you would also have variables which look
like:
\begin{boxverbatim}
TYPE(linked_list) :: head, tail
\end{boxverbatim}

Since Fortran pointers are not initialised in any particular state, you have
to remember to set the head and tail pointers to explicitly point nowhere
(normally called a null pointer by analogy with the older C style
pointers). This is done using the nullify command.
\begin{boxverbatim}
NULLIFY(head)
NULLIFY(tail)
\end{boxverbatim}

The same thing is important when creating a new linked list element, so you
would normally have a creation function for linked list elements.
\begin{boxverbatim}
SUBROUTINE create_element(element)

  TYPE(linked_list), POINTER :: element

  ALLOCATE(element)
  NULLIFY(element%next)
  NULLIFY(element%prev)

END SUBROUTINE create_element
\end{boxverbatim}
Note that the allocate function can be used on pointers in the same way that
it can be used with variables which have the allocatable attribute. However,
there is one important difference between a pointer and an allocatable
variable. If you attempt to allocate an already allocated variable which has
the allocatable attribute then the code will fail, whereas allocating an
already allocated pointer is perfectly valid, and will allocate the new
variable and point the pointer to it. This does not deallocate the memory that
the pointer previously pointed to, and Fortran does not have a ``garbage
collector'' which deallocates memory no longer accessible. So if you
allocate a pointer which already points to a variable, it is very important
that you have another pointer somewhere which points to the same memory. Once
you no longer have a pointer to an area of memory, that area of memory is
completely inaccessible and cannot even be deallocated. This is termed a
memory leak and for programs which run for many cycles and have a memory leak
on each cycle, the entire memory can very quickly be used up.

So, to add a new element to the list you would have a subroutine which looks
like:
\begin{boxverbatim}
SUBROUTINE add_element(element)

  TYPE(linked_list), POINTER :: element

  IF (.NOT. ASSOCIATED(head)) THEN
    ! Adding first element to list, so just set
    ! both head and tail to the element
    head=>element
    tail=>element
    RETURN
  ENDIF

  tail%next=>element
  element%prev=>tail
  tail=>element

END SUBROUTINE add_element
\end{boxverbatim}
This subroutine adds the new Fortran operator of ``$=>$'' which means ``points
to''. Unlike C or similar languages, Fortran pointers try to be partially
transparent to the end user, so the following code would fail:
\begin{boxverbatim}
PROGRAM test

  REAL, TARGET :: a = 10.0
  REAL, POINTER :: b

  b = a

END PROGRAM test
\end{boxverbatim}

This happens because Fortran will try to copy the value of ``a'' into ``b''.
However, ``b'' is a pointer which hasn't been initialised, so the code will
crash when it tries to copy the data in (in theory, the code may not crash if
the uninitialised ``b'' pointer happens to point somewhere in memory which is
a valid target, but this is very unlikely). Note also that ``a'' has the
attribute ``TARGET''. The target attribute means that it is possible to point
a pointer to this variable. You can only point a pointer to a variable which
is either a pointer itself or has the target attribute. This is to try and
keep Fortran pointers ``safer'' than C style pointers. The correct code would
use \inlinecode{b$=>$a}, at which point ``b'' is set to point to ``a'' and
can then be used everywhere in place of ``a''.

So, to set up a linked list of n elements, you would use the following code:
\begin{boxverbatim}
TYPE(linked_list), POINTER :: new
NULLIFY(new)

DO i = 1,n
  CALL create_element(new)
  CALL add_element(new)
ENDDO
\end{boxverbatim}

To then run through the elements of your newly created linked list, you would
use code like:
\begin{boxverbatim}
TYPE(linked_list), POINTER :: current

current=>head
DO WHILE(ASSOCIATED(current))
  ! Do stuff
  current=>current%next
ENDDO
\end{boxverbatim}

This code snippet introduces one new function ``ASSOCIATED'', which tells you
whether a pointer is a null pointer or not (this is why it is so important to
nullify new pointers, because ASSOCIATED on its own doesn't check whether a
pointer is valid, just whether or not it is a null pointer). You can also use
ASSOCIATED to check whether a pointer points to a particular object or not, in
which case the syntax is #RESULT = ASSOCIATED(b, TARGET=a)#, which
returns true if ``b'' points to ``a'', or false if it doesn't, even if ``b''
is a valid pointer pointing to something else. It also introduces the way in
which you must use linked lists in {\EPOCH}. The execution flow is as follows
\begin{itemize}
\item Point current to the current element to the start of the linked list
  (head).
\item Iterate while current points to a valid element.
\item Perform whatever actions you want on current.
\item Point current to the next element in the chain.
\end{itemize}
This leads to the slightly counter intuitive behaviour where even though the
loop only acts on the variable named ``current'', all of the elements in the
list are operated on. Although there are many tricks which can be performed
with linked lists, the only other aspect which needs to be explained is how
to delete elements. A subroutine to remove a single element from a linked list
would look like:
\begin{boxverbatim}
SUBROUTINE remove_element(element)

  TYPE(linked_list), POINTER :: element

  IF (ASSOCIATED(element%prev)) THEN
    ! Previous element exists
    element%prev%next=>element%next
  ELSE
    ! Previous element does not exist therefore element is the head. When
    ! element is removed the head is the element after the one being removed
    head=>element%next
  ENDIF

  IF (ASSOCIATED(element%next)) THEN
    ! next element exists
    element%next%prev=>element%prev
  ELSE
    ! next element does not exists therefore element is the tail. When element
    ! is removed the head is the element before the one being removed
    tail=>element%prev
  ENDIF

END SUBROUTINE remove_element
\end{boxverbatim}

Once again, this code looks slightly counter-intuitive, but if you go through
step by step, it's fairly simple. In the following discussion the element
being removed is called ``C'', the element before ``C'' (if it exists) is
called ``P'' and the element after ``C'' (if it exists) is called ``N''
\begin{itemize}
\item Check whether \inlinecode{C}'s prev element exists, this means that
  \inlinecode{P} exists.
\item If \inlinecode{P} exists then the element to be removed isn't at the
  start of the chain. When \inlinecode{C} is removed, we need
  \inlinecode{P}\%next to point to \inlinecode{C}\%next. This leads to the odd
  looking element\%prev\%next$=>$ element\%next syntax.
\item If \inlinecode{P} does not exist then \inlinecode{C} is at the the start
  of the chain. In order to not leave the chain orphaned when \inlinecode{C}
  is removed, we need head to point to \inlinecode{C}\%next.
\item Exactly the same logic applies for updating the element after
  \inlinecode{C}.
\item Check whether \inlinecode{C}'s next element exists, this means that
  \inlinecode{N} exists.
\item If \inlinecode{N} exists then the element to be removed isn't at the end
  of the chain. When \inlinecode{C} is removed, we need \inlinecode{N}\%prev
  to point to \inlinecode{C}\%prev.
\item If \inlinecode{P} does not exist then \inlinecode{C} is at the the start
  of the chain. In order to not leave the chain orphaned when \inlinecode{C}
  is removed, we need head to point to \inlinecode{C}\%next.
\end{itemize}

Therefore, code to remove some elements from a linked list would look like:
\begin{boxverbatim}
TYPE(linked_list), POINTER :: current, next

current=>head
DO WHILE(ASSOCIATED(current))
  next=>current%next
  IF (dealloc) THEN
    CALL remove_element(current)
    DEALLOCATE(current)
  ENDIF
  current=>next
ENDDO
\end{boxverbatim}
Note that ``current'' must be deallocated explicitly even after it has been
removed from the linked list to prevent a memory leak. Note also that the
pointer to the ``next'' element is saved before ``current'' is deallocated.
This is not necessary but means that there is only one IF statement rather
than the two otherwise required.

\subsubsection{Particles and particle species}
Particles are represented as linked lists of Fortran TYPES. The definition of
the particle type is:
\begin{boxverbatim}
  TYPE particle
    REAL(num), DIMENSION(3) :: part_p
    REAL(num), DIMENSION(c_ndims) :: part_pos
#ifdef PER_PARTICLE_WEIGHT
    REAL(num) :: weight
#endif
#ifdef PER_PARTICLE_CHARGE_MASS
    REAL(num) :: charge
    REAL(num) :: mass
#endif
    TYPE(particle), POINTER :: next, prev
#ifdef PARTICLE_DEBUG
    INTEGER :: processor
    INTEGER :: processor_at_t0
#endif
  END TYPE particle
\end{boxverbatim}
And the descriptions are
\begin{itemize}
\item #REAL(num) :: part_p(3)# - The particle momentum. Always dimension 3 even
  in 1D and 2D codes.
\item #REAL(num) :: part_pos(ndims)# - The particle position. Has
  the same dimensions as that of the code.
\item #REAL(num) :: weight# - The particle weight if the code is running with
  per particle weighting.
\item #REAL(num) :: charge# - The particle charge in Coulombs if the code is
  running with per particle charge.
\item #REAL(nun) :: mass# - The particle mass in kilograms if the code is
  running with per particle mass.
\item #TYPE(particle), POINTER :: next, prev# - The pointers to the next and
  previous elements of the linked list.
\item #INTEGER :: processor# - The rank of the processor that the particle
  thinks it is on. Used for debugging.
\item #INTEGER :: processor_at_t0# - The rank of the processor that the
  particle started on. Used for debugging.
\end{itemize}
Simply adding a new parameter to the definition of the particle type is NOT
sufficient to extend the particle type, since the communications when the
particle crosses a processor boundary do not know about the new parameter and
it will not be transmitted with the particle. How to add new properties to the
particle communication layer is described later.\\

The entire linked list of particles is encapsulated in another Fortran TYPE,
called \inlinecode{particle\_list}, which is defined as:
\begin{boxverbatim}
TYPE particle_list
  TYPE(particle), POINTER :: head
  TYPE(particle), POINTER :: tail
  INTEGER(KIND=8) :: count
  ! Pointer is safe if the particles in it are all unambiguously linked
  LOGICAL :: safe

  TYPE(particle_list), POINTER :: next, prev
END TYPE particle_list
\end{boxverbatim}
And its properties are:
\begin{itemize}
\item #TYPE(particle), POINTER :: head# - The first particle in the linked list.
\item #TYPE(particle), POINTER :: tail# - The last particle in the linked
  list. New particles added to the end of the list are added onto the end of
  the tail element, and the new last particle becomes the new tail element.
\setlength{\emergencystretch}{3em}
\item #INTEGER(KIND=8) :: count# - The number of particles in this particle
  list. Note that the \inlinecode{particle\_list} type is not
  directly MPI aware, so this is literally the number of particles in
  {\it this} particle list, not the number of particles of this species on all
  processors.
\item #LOGICAL :: safe# - A particle list is {\it safe} if the particles in it
  are unambiguously linked. That is that the \inlinecode{count}th particle is
  guaranteed to have its \inlinecode{next} property be null. Most particle
  lists within {\EPOCH} are safe, but sometimes it is useful to be able to have
  particle lists which are subsets of longer particles lists, and these
  particle lists are not {\it safe}.
\item #TYPE(particle_list), POINTER :: next, prev# - At present, {\EPOCH} does
  not use these pointers, which are intended to allow multiple particle lists to
  be attached together. Certain parts of {\EPOCH}, such as the I/O system are
  aware of these pointers and will automatically use them if they are ever
  set. They are reserved for future use.
\end{itemize}

The \inlinecode{particle\_list} objects are used to abstract all the functions
of the linked list, including adding and removing particles and transporting
particles between processors.\\
\\
The particle species are represented by yet another Fortran TYPE, this time
called \inlinecode{particle\_species}, which is defined as:
\begin{boxverbatim}
  TYPE particle_species
    ! Core properties
    CHARACTER(string_length) :: name
    TYPE(particle_species), POINTER :: next, prev
    INTEGER :: id
    INTEGER :: dumpmask

    REAL(num) :: charge
    REAL(num) :: mass
    REAL(num) :: weight
    INTEGER(KIND=8) :: count
    TYPE(particle_list) :: attached_list

#ifdef TRACER_PARTICLES
    LOGICAL :: tracer
#endif

    ! particle cell division
#ifdef SPLIT_PARTICLES_AFTER_PUSH
    INTEGER(KIND=8) :: global_count
    LOGICAL :: split
    INTEGER(KIND=8) :: npart_max
    ! Secondary list
    TYPE(particle_list), DIMENSION(:,:), POINTER :: secondary_list
#endif

    ! Injection of particles
    INTEGER(KIND=8) :: npart_per_cell
    REAL(num), DIMENSION(:), POINTER :: density
    REAL(num), DIMENSION(:,:), POINTER :: temperature

    ! Species_ionisation
#ifdef PARTICLE_IONISE
    LOGICAL :: ionise
    INTEGER :: ionise_to_species
    INTEGER :: release_species
    REAL(num) :: critical_field
    REAL(num) :: ionisation_energy
#endif
    ! Attached probes for this species
#ifdef PARTICLE_PROBES
    TYPE(particle_probe), POINTER :: attached_probes
#endif
  END TYPE particle_species
\end{boxverbatim}
Again, most of these properties are self explanatory, but they are detailed
below.
\begin{itemize}
\item #CHARACTER(LEN=entry_length) :: name# - The name of the particle
  species. Used when constructing things like ``ekbar\_electron'' and similar
  names.
\item #TYPE(particle_species), POINTER :: next, prev# - Particle species are
  connected to each other as a linked list using pointers as well as being
  available through a simple array. These pointers are used behind the scenes
  in the I/O.
\item #INTEGER :: id# - The number of the species, so for the species
  \inlinecode{species\_list(1)}, the id field would be 1. For
  \inlinecode{species\_list(2)}, the id field would be 2 etc.
\item #INTEGER :: dumpmask# - Bitmask to determine when this species should be
  dumped in diagnostic output.
\item #REAL(num) :: charge# - The charge on a single particle of the species in
  Coulombs.
\item #REAL(num) :: mass# - The mass of a single particle of the species in
  kilograms.
\item #REAL(num) :: weight# - The per-species particle weight.
\item #INTEGER(KIND=8) :: count# - The number of particles of this species on
  all processors. NOTE that this is only accurate if the code is compiled with
  the correct preprocessor options. Without the correct preprocessor options,
  this will be accurate at the start of the code runtime, but will not be if
  any particles enter or leave the domain. This is mainly a debugging
  parameter.
\setlength{\emergencystretch}{3em}
\item #TYPE(particle_list) :: attached_list# - This is the
  \inlinecode{particle\_list} object which holds the particles assigned to this
  species on this processor. Particles are attached to this list at all
  times EXCEPT when they are explicitly split when the code is compiled with
  the \inlinecode{SPLIT\_PARTICLES\_AFTER\_PUSH} option. When the
  code is compiled with this option, the particles are attached to
  \inlinecode{attached\_list} except between the calls to
  \inlinecode{reorder\_particles\_to\_grid} and
  \inlinecode{reattach\_particles\_to\_mainlist} in
  \inlinecode{epoch{\it n}d.F90} where the particles are instead attached to
  \inlinecode{secondary\_list}. This is explained later.
\item #LOGICAL :: tracer# - Whether or not this species is a tracer particle. If
  a species is a tracer species then it moves under the fields as normal for a
  particle with its mass and charge but contributes no current.
\item #LOGICAL :: split# - {\EPOCH} includes a very early version of a particle
  splitting operator. It works mechanically but has undesirable properties at
  present. If this flag is true then the code attempts to split the particles
  when the pseudoparticle number density drops too low.
\item #INTEGER(KIND=8) :: npart_max# - Used with the particle splitting
  operator. When the total number of particles equals this number, further
  particle splitting is suppressed.
\item #TYPE(particle_list), DIMENSION(:,:,:), POINTER :: secondary_list# - When
  the code is compiled with the \inlinecode{-DSPLIT\_PARTICLES\_AFTER\_PUSH}
  option, the code allocates
  \inlinecode{secondary\_list(0:nx+1,0:ny+1,0:nz+1)} and then loops over all
  particles. It calculates the cell in which each particle is and moves the
  particle from \inlinecode{attached\_list} to the correct element of
  \inlinecode{secondary\_list} for that cell. This means the particles which
  are nearby in space are now linked together in an array of linked lists.
  This allows things such as collision operators which require direct
  interaction between nearby particles.
\item #INTEGER(KIND=8) :: npart_per_cell# - The number of pseudoparticles per
  cell in the initial conditions. This is used with the moving window function
  to ensure that the same number of particles per cell are used for the new
  material introduced at the leading edge of the window.
\item #REAL(num), DIMENSION(:,:), POINTER :: density# - The density of the
  plasma at the leading edge of the window at the start of the simulation. This
  is used to structure the density of the new material introduced at the leading
  edge of the plasma.
\item #REAL(num), DIMENSION(:,:,:), POINTER :: temperature# - The temperature
  in of the plasma at the leading edge of a moving window at the start of the
  simulation. The final index of the array is the direction in which the
  temperature is set (1=x, 2=y, 3=z).
\item #LOGICAL :: ionise# - If the ionisation model is activated then this
  species should ionise.
\item #INTEGER :: ionise_to_species# - The species number for the next ionised
  state of this species.
\item #INTEGER :: release_species# - Specifies what type of particle should be
  released when this species ionises (i.e. which species is the electron).
\item #REAL(num) :: ionisation_energy# - The ionisation energy for the next
  ionisation of this species.
\item #TYPE(particle_probe), POINTER :: attached_probes# - A pointer pointing to
  the head of an attached linked list of particle probe diagnostics.
\end{itemize}

\subsubsection{EM Fields}
There are nine variables which are used in updating the EM field solver. These
are
\begin{itemize}
\item ex - Electric field in the X direction.
\item ey - Electric field in the Y direction.
\item ez - Electric field in the Z direction.
\item bx - Magnetic field in the X direction.
\item by - Magnetic field in the Y direction.
\item bz - Magnetic field in the Z direction.
\item jx - Current in the X direction.
\item jy - Current in the Y direction.
\item jz - Current in the Z direction.
\end{itemize}
The EM fields in {\EPOCH} are simple allocatable arrays, which are of size
(-2:nx+3,-2:ny+3,-2:nz+3), although this includes the ghost cells. The length of
the core domain is different for each variable due to the grid stagger.

The {\EPOCH} field solver is a Yee staggered 2nd order FDTD scheme, directly
based on the scheme in the PSC by Hartmut Ruhl and is contained in the file
\inlinecode{fields.f90}. To locate a variable on the grid there is a simple
rule.
\begin{itemize}
\item Start at the cell centre.
\item For an $E$ field component, move the field half a grid point in the
  direction that the field points if possible.
\item For a $B$ field component, move the field half a grid point in all
  directions {\it except} the one it points.
\end{itemize}
This is illustrated in Figure~\ref{yeegrid} for the 2D case.\\

\captionedimage{./images/stagger}{yeegrid}{The Yee grid in 2D}


The grid stagger means that you have to be careful with boundary conditions
since some variables are defined on the domain boundaries whereas others are
defined on either side of a domain boundary. This is handled automatically by
the built in boundary routines, but must be understood if developing other
boundary conditions. To explain it, consider only the left/right boundary in 1D
and consider $E_x$ and $B_x$.\\

$E_x$ is defined on the cell boundary, so \inlinecode{ex(0)} is the value of
$E_x$ on the left boundary and similarly \inlinecode{ex(nx)} is the
value on the right boundary. Conversely, in the 1D code $B_x$ is cell centred
(in reality, $B_x$ is never used in the field update and is unimportant since
any gradients in $B_x$ in 1D automatically break the solenoidal condition, but
this is still a useful example.). This means that \inlinecode{bx(1)} is the
centre of the first cell in the domain, and \inlinecode{bx(0)} is the value at
the centre of the first left hand ghost cell. This means the you must do
different things as boundary conditions for the two fields for some boundary
conditions.\\

For example, if you want to clamp the value of $E_x$ to be zero on the
boundary, then just set \inlinecode{ex(0) = 0.0\_num} since \inlinecode{ex(0)}
lies on the boundary. To do the same for $B_x$ on the boundary you have to
set \inlinecode{bx(0) = -bx(1)}. This is because if you use a linear
reconstruction of $B_x$ (i.e second order) then the point between
\inlinecode{bx(0)} and \inlinecode{bx(1)} has the value
$B_x(1/2) = \left(B_x(1)+B_x(0)\right)/2$. Similarly, if you want to set zero
gradient on the boundary then for $E_x$ you set \inlinecode{ex(-1) = ex(1)},
whereas for $B_x$ you would set \inlinecode{bx(0) = bx(1)}. This is explained
in more detail in \sect{bcs}.

In the particle pusher, time centred field variables are needed for second
order accuracy, so an FDTD scheme is used to advance the fields. This looks
like

\begin{itemize}
\item $\vec{E}^{n+\frac{1}{2}} = \vec{E}^n + \frac{\Delta t}{2} \left( c^2
  \nabla \wedge \vec{B}^{n} -\vec{j}^{n} \right)$
\item $\vec{B}^{n+\frac{1}{2}} = \vec{B}^n - \frac{\Delta t}{2} \left( \nabla
  \wedge \vec{E}^{n+\frac{1}{2}} \right)$
\item Call particle pusher which calculates $j^{n+1}$ currents
\item $\vec{B}^{n+1} = \vec{B}^{n+\frac{1}{2}} - \frac{\Delta t}{2} \left(
  \nabla \wedge \vec{E}^{n+\frac{1}{2}} \right)$
\item $\vec{E}^{n+1} = \vec{E}^{n+\frac{1}{2}} + \frac{\Delta t}{2} \left( c^2
  \nabla \wedge \vec{B} ^{n+1} - \vec{j}^{n+1} \right)$
\end{itemize}
Note that all spatial derivatives are calculated using the staggered grid, so
the final derivatives in the code appear one sided. However, this is not the
case, and all spatial derivatives are second order accurate. Higher order
spatial derivatives schemes for {\EPOCH} are being developed to improve the
dispersion properties of the code when resolving small timescales.

\subsection{The particle pusher}
{\EPOCH}'s particle pusher is based on the one from the PSC by Hartmut Ruhl, and
is a Birdsall and Landon type PIC scheme using Villasenor and Buneman current
weighting. It is contained in the file \inlinecode{particles.F90}. The
operation of the particle pusher is fairly simple, but there are a few elements
which need some clarification.
\begin{itemize}
\item The update to the particle momenta etc. does not explicitly include the
  particle weight function. This means that the pseudoparticle momenta etc. are
  the momentum for a single real particle of the collection of real particles
  represented by that pseudoparticle, NOT the momentum of the whole collection
  of real particles.
\item \inlinecode{gamma} - The variable gamma which appears in various places
  is the relativistic $\gamma$ which is needed to convert the particle momentum
  into the particle velocity using the relation $\vec{p} = \gamma m \vec{v}$.
  Here, $\vec{p}$ is the particle momentum, $\vec{v}$ is the particle
  velocity and $m$ is the particle rest mass.
  If {\EPOCH} was not relativistic then this would simply be $1$. Since
  {\EPOCH} is relativistic, gamma is defined as
  $\left(\vec{p}.\vec{p}/m c^2 + 1\right)^{1/2}$.
\item \inlinecode{cell\_x1=cell\_x1+1} - There are lines like this after all
  the sections of the routine where the cell a particle is in is
  calculated. This is because, for a cell centred variable, the domain runs
  (1:nx,1:ny,1:nz) rather than (0:nx-1,0:ny-1,0:nz-1).
\end{itemize}

\subsubsection{Particle shape functions}
The key feature of a PIC code controlling the smoothness of the solution is the
particle shape function. That is the function that describes the assumed
distribution of the real particles making up a pseudoparticle. The simplest
solution is to assume that the pseudoparticles uniformly fill the cell in which
the pseudoparticle is located. This has the advantages of speed and simplicity
but produces very noisy solutions. The next simplest approach is to assume a
triangular shape function with the peak of the triangle located at the position
of the pseudoparticle and a width of $ 2 \Delta x$, as illustrated in
Figure~\ref{shape}. This is the approach used
in {\EPOCH} and is a good trade-off between cleanness of solution and
speed. Higher order methods based on spline interpolation can be used and do
produce smoother solutions, but they are significantly slower and the benefits
of the schemes can easily be overstated. {\EPOCH} does now include an option to
use 3rd order b-spline interpolation in all parts of the code. This option is
enabled with the \inlinecode{-DPARTICLE\_SHAPE\_BSPLINE3} compile time option
in the makefile.\\
\[
S(w) =
\begin{cases}
1 - \frac{x_i-w}{\Delta x}, & (x_i-w) \le \Delta x \\
0, & \mbox{otherwise}
\end{cases}
\]

Functions derived from the particle shape function appear in two places in the
core solver: when the EM fields are interpolated to the position of the
pseudoparticle and when the current is updated and properties of the
pseudoparticle are copied onto the grid. These two uses of the shape function
are conceptually similar, but have different forms.

\captionedimage{./images/shape}{shape}{Second order particle shape function}

To derive the equations for calculating the field acting on a particle,
you calculate the overlap of the particle shape function with the function
representing the fields on the grid. In {\EPOCH}, the fields are approximated at
first order so that the field is constant over each cell. Consider a particle
with position $X$, where $X$ lies in the cell centred at $x_i$ and grid
spacing $\Delta x$. The integral is split into four parts; that part of the
shape function which overlaps with the cell $x_{i-1}$, the part of the shape
function from the left boundary of $x_i$ to the point of the triangle, the part
of the shape function from the point of the triangle to the right hand edge of
$x_i$ and finally that part of the shape function which overlaps cell
$x_{i+1}$. Assuming that fields are constant inside each cell, this takes the
form
\[
\begin{aligned}
  Fpart_i~=~\frac{1}{\Delta x} & \left[ \int^{x_{i-1}
      + \frac{\Delta x}{2}}_{X-\Delta x}
      F_{i-1} \left( 1-\frac{X-x}{\Delta x} \right) dx \right.\\
    & + \int^{X}_{x_i-\frac{\Delta x}{2}}
      F_i \left( 1-\frac{X-x}{\Delta x} \right) dx \\
    & + \int^{x_i+\frac{\Delta x}{2}}_{X}
      F_i \left(  1-\frac{x-X}{\Delta x}\right) dx \\
    & \left. + \int^{X+\Delta x}_{x_{i+1} - \frac{\Delta x}{2}}
      F_{i+1} \left( 1-\frac{x-X}{\Delta x}\right) dx\right]
\end{aligned}
\]

Performing these integrals and remembering that $x_{i-1}+\frac{\Delta x}{2}$ is
equal to $x_i-\frac{\Delta x}{2}$ since the grid is uniformly spaced with
spacing $\Delta x$ this gives a final formula for the field at a particle of

\[
\begin{aligned}
  F_{part}~=~& \frac{1}{2} F_{i-1} \left( \frac{1}{2}
  + \frac{x_i-X}{\Delta x} \right)^2 \\
  & + F_i \left( \frac{3}{4} - \frac{(x_i-X)^2}{\Delta x^2} \right)\\
  & +\frac{1}{2} F_{i+1} \left( \frac{1}{2} - \frac{x_i-X}{\Delta x} \right)^2
\end{aligned}
\]

%If you are running the code with the \inlinecode{-DSPLINE\_FOUR} high order
%particle weight function option then the triangular shape function is replaced
%by a 4th order spline function which has a basis of 5 cells rather than 3
%cells. The form of this function is

In the code calculating the strength of a cell centred field on the particle
is done as follows.
\begin{boxverbatim}
  REAL(num) :: cell_x_r, cell_frac_x
  INTEGER :: cell_x1
  REAL(num) :: gx(-2:2)
  TYPE(particle), POINTER :: current

  part_x  = current%part_pos(1) - x_min_local

  ! Work out the grid cell number for the particle.
  ! Not an integer in general.
  cell_x_r = part_x / dx

  ! Round cell position to nearest cell
  cell_x1 = FLOOR(cell_x_r + 0.5_num)
  ! Calculate fraction of cell between nearest cell boundary and particle
  cell_frac_x = REAL(cell_x1, num) - cell_x_r
  cell_x1 = cell_x1 + 1

  cf2 = cell_frac_x**2
  gx(-1) = 0.25_num + cf2 + cell_frac_x
  gx( 0) = 1.5_num - 2.0_num * cf2
  gx( 1) = 0.25_num + cf2 - cell_frac_x

  f_part = &
        gx(-1) * F(cell_x1-1) &
      + gx( 0) * F(cell_x1  ) &
      + gx( 1) * F(cell_x1+1)
\end{boxverbatim}

where \inlinecode{f\_part} is the field at the particle location. Note that
this has been simplified a little for brevity. Just the triangle shape
function is given.
In 2D or 3D, you just calculate gy
in the same manner as gx and calculate the weight over all the cells affected
by the individual 1D shape functions. In 2D this looks like:
\begin{boxverbatim}
  f_part = 0.0_num
  DO iy = sf_min, sf_max
    DO ix = sf_min, sf_max
      f_part = f_part + f(cell_x+ix, cell_y+iy) * gx(ix) * gy(iy)
    ENDDO
  ENDDO
\end{boxverbatim}
The variables \inlinecode{sf\_min} and \inlinecode{sf\_max} contain the shape
function order parameters which indicate the cells each side of the cell
containing the particle which are overlapped by the particle shape function.
They are defined in \inlinecode{shared\_data.F90} and should only be changed
by the developer if a new particle shape function is being added.
Although provided here as pseudo-code for the particle push, it should be
noted that the actual particle push unrolls these loops for the sake of speed.

Inside the particle pusher the $E$ and $B$ fields are not cell centred fields,
but Yee staggered. This means that there is a small change to the above
mentioned example. In 1D this change looks like

\begin{boxverbatim}
  REAL(num) :: cell_x_r, cell_frac_x
  INTEGER :: cell_x1, cell_x2
  REAL(num) :: gx(-2:2), hx(-2:2)
  TYPE(particle), POINTER :: current

  part_x  = current%part_pos(1) - x_min_local

  ! Work out the grid cell number for the particle.
  ! Not an integer in general.
  cell_x_r = part_x / dx

  ! Round cell position to nearest cell
  cell_x1 = FLOOR(cell_x_r + 0.5_num)
  ! Calculate fraction of cell between nearest cell boundary and particle
  cell_frac_x = REAL(cell_x1, num) - cell_x_r
  cell_x1 = cell_x1 + 1

  ! Calculate weights
  INCLUDE 'include/triangle/gx.inc'

  ! Now redo shifted by half a cell due to grid stagger.
  ! Use shifted version for ex in X, ey in Y, ez in Z
  ! And in Y&Z for bx, X&Z for by, X&Y for bz
  cell_x2 = FLOOR(cell_x_r)
  cell_frac_x = REAL(cell_x2, num) - cell_x_r + 0.5_num
  cell_x2 = cell_x2 + 1

  ! Calculate weights
  INCLUDE 'include/triangle/hx_dcell.inc'

  ! bx is cell centred
  bx_part = 0.0_num
  DO ix = sf_min, sf_max
    bx_part = bx_part + bx(cell_x1+ix) * gx(ix)
  ENDDO

  ! ex is staggered 1/2 a cell to the right
  ex_part = 0.0_num
  DO ix = sf_min, sf_max
    ex_part = ex_part + ex(cell_x2+ix) * hx(ix)
  ENDDO
\end{boxverbatim}

In 2D and 3D, you just combine the shifted and unshifted shape functions and
associated cell positions depending on the position of the variable in the
cell. Therefore, in 3D and using the loop notation for clarity you would get:
\begin{boxverbatim}
  DO iz = sf_min, sf_max
    DO iy = sf_min, sf_max
      DO ix = sf_min, sf_max
        ex_part = ex_part + hx(ix) * gy(iy) * gz(iz) * &
            ex(cell_x2+ix, cell_y1+iy, cell_z1+iz)
      ENDDO
    ENDDO
  ENDDO

  DO iz = sf_min, sf_max
    DO iy = sf_min, sf_max
      DO ix = sf_min, sf_max
        bx_part = bx_part + gx(ix) * hy(iy) * hz(iz) * &
            bx(cell_x2+ix, cell_y1+iy, cell_z1+iz)
      ENDDO
    ENDDO
  ENDDO
\end{boxverbatim}

Since $E_x$ is staggered half a grid cell in the x direction, whereas $B_x$ is
staggered by half a grid cell in the y and z directions.

The next stage is to consider how to copy pseudoparticle
properties on the grid. This is very similar to the function for calculating
grid variables at the particle location and, for each grid point $x_i$,
consists of integrating the part of the particle shape function which overlaps
the $i^{th}$ cell. That is

\[
F(i) = Data  \int^{x_i+\frac{\Delta x}{2}}_{x_i-\frac{\Delta x}{2}} S(X-x) dx
\]

Where $Data$ is the particle property to be copied onto the grid. In {\EPOCH},
since the particle shape function is known to go to zero outside a distance of
$2 \Delta x$ from the maximum, the maximum number of cells that can possibly be
overlapped by a given particle shape function is 3; the cell containing the
particle maximum and the two cells to either side. Performing the integration
using the triangular shape function given above gives the result

\[
  F(i) =
\begin{cases}
  \frac{3}{4} - \frac{|X-x_i|^2}{\Delta x^2}, & |X-x_i|
    \le \frac{\Delta x}{2}\\
  \frac{1}{2} \left(\frac{3}{2}
    - \frac{|X - x_i|}{\Delta x} \right)^2, & \frac{\Delta x}{2} < |X-x_i|
    \le \frac{3 \Delta x}{2}\\
  0, & |X-x_i| > \frac{3 \Delta x}{2}\\
\end{cases}
\]

%Again, when using high order particle shape functions using the
%\inlinecode{-DSPLINE\_FOUR} option this is replaced with an equivalent form for
%a 4th order spline.

When this is translated into the code, it looks very similar to that presented
for the case where grid properties are interpolated to the particle
position. This form is used in the particle pusher to perform the current
update and in the routines in \inlinecode{src/io/calc\_df.F90} to copy particle
properties onto the grid for output. The form from calc\_df is rather clearer
and easier to see in operation. In 1D it looks like:
\begin{boxverbatim}
  cell_x_r = (current%part_pos - x_min_local) / dx + 1.5_num
  cell_x = FLOOR(cell_x_r)
  cell_frac_x = REAL(cell_x, num) - cell_x_r + 0.5_num

  CALL particle_to_grid(cell_frac_x, gx)

  wdata = part_m * fac
  DO ix = sf_min, sf_max
    data_array(cell_x+ix) = data_array(cell_x+ix) + gx(ix) * wdata
  ENDDO
\end{boxverbatim}

Once again multi-dimensional codes just have the weighting functions multiplied
together.
\begin{boxverbatim}
  DO iy = sf_min, sf_max
    DO ix = sf_min, sf_max
      data_array(cell_x+ix, cell_y+iy) = &
          data_array(cell_x+ix, cell_y+iy) + gx(ix) * gy(iy) * wdata
    ENDDO
  ENDDO
\end{boxverbatim}

\subsubsection{Current calculation}
{\EPOCH} uses the Villasenor and Buneman (Villasenor and Buneman, Computer
Physics Communications 69(1992) 306-316) current calculating scheme which
solves the additional equation
${\partial \rho}/{\partial t} = \nabla\cdot\vec{J}$ to
calculate the current at each timestep. The main advantage of this scheme is
that it conserves charge {\it on the grid} rather than just globally conserving
charge on the particles. This means that the error in the solution to Gauss's
law is conserved, so if Gauss's law is satisfied for $t = 0$ it
remains satisfied for all time.\\

The Villasenor and Buneman scheme works because exactly the same charge added
to one cell is subtracted from another cell, which in turn means that exactly
the same current added to one cell is subtracted from another cell. This is
intuitively correct since a point particle crossing a cell boundary would
represent the loss of that particle's contribution to the current from the
source cell and the gain of that particle's contribution to the current by the
destination cell. In fact this simple type of cell boundary crossing
current calculation was used in classical Buneman type PIC codes.\\

The scheme is messy, in practise, but simple. After the main particle push, the
particle is advanced a further half timestep into the future to first order
using the velocities calculated at the end of the particle push. The particle
position at $t + dt/2$ were stored earlier, and combined with the newly
calculated particle position at $t + {3dt}/{2}$ this allows a time centred
evaluation of ${\partial \rho}/{\partial t}$ meaning that the current
update is second order accurate in time. The spatial order of the scheme
matches the spatial order of the particle weight function.\\

The weight functions for transferring particle properties onto the grid at the
two timesteps are calculated including a shift when necessary to allow for the
particle having crossed a cell boundary. Since the charge associated with the
particle is spatially distributed using the weight function, all that is
necessary to calculate ${\partial \rho}/{\partial t}$ is to subtract the
two functions, multiply by the charge on the pseudoparticle and the
pseudoparticle weight and finally divide by $dt$. The spatial derivative of
$\vec{J}$ is then converted to a one sided finite difference form and solved
directly. In multiple dimensions this is slightly complicated by the effects of
offsets in directions other than the direction that a given current component
is pointing in, with this adding additional weight factors based on the overlap
of the shape functions in other directions. This is explained in full in the
Villasenor and Buneman paper already quoted.\\

Currents in ignorable directions are simply calculated using $J = n\rho\vec{v}$
with the correct shape functions to ensure that the current is placed in the
correct places.

\subsection{The timestep}

The timestep is calculated in the subroutine \inlinecode{set\_dt} in the file
\inlinecode{src/housekeeping/setup.F90}. All that the subroutine has to do is set
the variable \inlinecode{dt} to set the timestep for the whole code. Any
additional timestep constraints should be coded into this subroutine. This
should be implemented after the existing \inlinecode{dt=} lines but before the
line \inlinecode{dt = dt\_multiplier * dt}. Such a modification should be set so
that it only changes the timestep if the timestep is MORE restrictive than that
calculated from the core code. An example would be:
\begin{boxverbatim}
  dt = dx * dy / SQRT(dx**2 + dy**2) / c
  dt = MIN(dt, my_new_dt)
\end{boxverbatim}

In the core {\EPOCH} code the timestep can be calculated identically on each
processor, so there is no requirement to synchronise the timestep across
multiple processors. If your new timestep restriction uses information local to
each processor then some additional lines must be added to the
\inlinecode{set\_dt} routine after the timestep has been calculated which
should read:
\begin{boxverbatim}
  REAL(num) :: dt_global
          .
          .
          .
  CALL MPI_ALLREDUCE(dt_global, dt, 1, mpireal, MPI_MIN, comm, errcode)
  dt = dt_global
\end{boxverbatim}

This uses another MPI command to determine the most restrictive timestep across
all processors. {\EPOCH} is not written in a way that permits operation with
different timesteps on different processors, and the behaviour of the code is
undefined (and likely wrong) if the code runs with different timesteps on
different processors.

\subsection{Boundary conditions}
\label{sec:bcs}
Boundary conditions in {\EPOCH} are split into three types
\begin{itemize}
\item Simple field boundaries.
\item Laser and outflow boundaries.
\item Particle boundaries.
\end{itemize}

These boundaries can be combined in different ways to give different
effects. From the end user perspective there are 4 boundaries which can be
applied to each edge of the simulation domain. These are
\begin{itemize}
\item Periodic
  \subitem Particles periodic
  \subitem Fields periodic
  \subitem Lasers off
\item Other
  \subitem Particles reflect
  \subitem Fields clamped zero
  \subitem Lasers off
\item Simple Laser
  \subitem Particles transmissive
  \subitem Fields clamped zero
  \subitem Lasers applied at half timestep for $B$ field
\item Simple outflow
  \subitem Particles transmissive
  \subitem Fields clamped zero
  \subitem No lasers applied at half timestep, but outflow conditions applied
  to $B$ field at half timestep
\end{itemize}
The boundary conditions are applied in too many places in the code to give a
full description of them, but the laser boundaries are only applied in
\inlinecode{src/fields.f90}. The boundaries requested by the user are converted
into the conditions on the fields and particles in the routine
\inlinecode{setup\_particle\_boundaries} in
\inlinecode{src/boundaries.F90}. For each of the six possible boundaries
(x\_min, x\_max, y\_min, y\_max, z\_min, z\_max) there is a variable which will
be named something like \inlinecode{bc\_x\_min\_particle} or
\inlinecode{bc\_y\_max\_field} which controls the boundary condition which will
be applied to either the field or the particles.

\subsubsection{Simple field boundaries}
There are two subroutines which apply the standard boundary conditions:
\inlinecode{field\_zero\_gradient} and \inlinecode{field\_clamp\_zero}. The
type of boundary condition that the two apply is obvious from the name, but
the two functions have different calling conventions.

\pagebreak
\begin{codedef}
SUBROUTINE field_zero_gradient
\HRule
REAL(num), DIMENSION(-2:,-2:,-2:), INTENT(INOUT) :: field
INTEGER, INTENT(IN) :: stagger_type, boundary
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/boundary.f90}\\
\\[0.5cm]
{\Large Description\\}
#field_zero_gradient# is the routine which applies zero gradient
boundary conditions to a field variable passed in the parameter
\inlinecode{field}. The remaining two parameters are the stagger type and
the boundary number. These are named constants defined in
\inlinecode{src/shared\_data.F90}, \inlinecode{c\_stagger\_*} for the
stagger type and \inlinecode{c\_bd\_*} for the boundary number.
The routine can be used as a global field boundary condition by
setting one of the field boundary conditions to c\_bc\_zero\_gradient in
\inlinecode{setup\_particle\_boundaries}, but it is mostly used to give
boundary conditions for the autoloader.
\\[0.5cm]
{\Large Notes\\ \\}

\pagebreak
\begin{codedef}
SUBROUTINE field_clamp_zero
\HRule
REAL(num), DIMENSION(-2:,-2:,-2:), INTENT(INOUT) :: field
INTEGER, INTENT(IN) :: stagger_type, boundary
\end{codedef}
\vspace{1cm}
{\Large Description\\}
\inlinecode{field\_clamp\_zero} is the routine which clamps the field given by
the \inlinecode{field} variable to zero on the boundary.
The remaining two parameters are the stagger type and
the boundary number. These are named constants defined in
\inlinecode{src/shared\_data.F90}, \inlinecode{c\_stagger\_*} for the
stagger type and \inlinecode{c\_bd\_*} for the boundary number.
\\[0.5cm]
{\Large Notes\\}
\pagebreak

Additional boundary conditions should follow the same basic principle as these
routines. Note that all of the routines test for \inlinecode{\{x,y,z\}\_\{min,max\}\_boundary} to confirm that a given processor is at the
edge of the real domain, and so should have a real boundary condition applied
to it. This also explains why there are no explicit periodic boundary condition
routines, since by connecting the processors cyclically in a periodic direction
the domain boundary effectively becomes another internal processor boundary.

\subsubsection{Laser and outflow boundaries}
Setting up lasers is explained in \sect{laser_block}
and the way in which the laser specification in the input deck works is
described later, so this section just describes how the actual laser and
outflow routines in \inlinecode{laser.f90} work.\\

The laser boundaries in {\EPOCH} are based on a rewriting of Maxwell's
equations in a new form which expresses the fields explicitly in terms of waves
propagating in both directions along each co-ordinate axis with both S and P
polarisation states. In the X direction,
\[
\partial_t(E_y \pm B_z) \pm \partial_x(E_y \pm B_z) = \pm \partial_yE_x
+ \partial_zB_x -\frac{j_y}{\epsilon_0}
\]
\[
\partial_t(E_z \mp B_y) \pm \partial_x(E_z \mp B_y) = \pm \partial_zE_x
- \partial_yB_x -\frac{j_z}{\epsilon_0}
\]
It is then possible to rewrite these equations to provide a boundary condition
on $B_z$ and $B_y$ to give propagating EM waves at the boundary. For waves
travelling into the boundary, this gives a transmissive boundary, and if the
components for waves propagating out from the boundary are set to be non-zero
then it also introduces an EM wave propagating from the left boundary.\\

This boundary condition is found in the file \inlinecode{laser.f90} which also
includes the routines for handling the \inlinecode{laser\_block} objects which
represent how lasers are represented in {\EPOCH}.

\subsubsection{Particle boundaries}
Due to the time that is required to loop over all the particles the particle
boundary conditions in {\EPOCH} combine the inter-processor boundary conditions
with the real boundary conditions. The boundary conditions for particles are in
the routine \inlinecode{particle\_bcs} in the file \inlinecode{boundary.f90} \\
Currently {\EPOCH} includes only three particle boundary conditions
\begin{itemize}
\item c\_bc\_open - Particles pass through the boundary and are destroyed. Total
  pseudoparticle number is not conserved in this mode.
\item c\_bc\_periodic - Particles which leave one side of the box reappear on
  the other side.
\item c\_bc\_reflect - Particles reflect off the boundary as if it was a hard
  boundary.
\end{itemize}

Although the routine looks rather messy, it is fairly easy to understand. The
sequence goes:
\begin{itemize}
\item Loop over all species.
\item Create particle list objects for particles to be sent to and received from
  other processors.
\item Loop over all particles in the species.
\item If the particle has crossed a local boundary then it either
  has crossed the boundary of the problem domain and needs a boundary condition
  to be applied or it has just crossed a processor boundary and
  needs to be transferred to a neighbouring process. Set
  \inlinecode{\{x,y,z\}bd} which is used to identify which processor relative
  to the current processor the particle potentially needs to be moved to and
  then test for the domain boundary.
\item If the particle has crossed an open domain boundary then either add it to
  another list to be dumped to disk if the user has requested this, or
  otherwise just deallocate the particle to reclaim memory. Set
  \inlinecode{out\_of\_bounds} to indicate that the particle has left the
  system.
\item If the particle has crossed the domain boundary and that boundary has
  reflecting boundary conditions then reflect the particle.
\item If the particle has crossed the domain boundary and that boundary has
  periodic boundary conditions then move the particle to the opposite side
  of the domain.
\item End particle loop.
\item Remove all the particles which have left the current process.
\item Loop over all possible neighbouring processors for the current processor
  and exchange particle lists with that processor.
\item Add any received particles onto the particle list for the current
  species.
\item End species loop.
\end{itemize}

Note that, unlike for fields, there is explicit periodic boundary code. This is
because although the MPI routines place the particle on the correct processor
after the MPI routines, the particle's position variable still places it beyond
the other end of the domain. The MPI parallelism for exchanging particles is
hidden in the routines which deal with the particle list objects and are
described in the next section.
\pagebreak

\subsubsection{MPI Boundaries}
There are three routines which deal with MPI exchange for field variables in
{\EPOCH}. Two are closely related and will be considered together. The third
deals with using MPI to sum variables at processor boundaries rather than
synchronise ghost cells.\\\\

\begin{codedef}
SUBROUTINE field_bc
\HRule
REAL(num), DIMENSION(-2:,-2:,-2:), INTENT(INOUT) :: field
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/boundary.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{field\_bc} exchanges the information in the ghost cells between
adjacent processors. Any field variable which is used in a calculation that
requires operations involving information from points other than the
current point should call this routine each time the variable is updated. This
will ensure that the ghost cells are populated from adjacent processors.
(i.e. if you only need to access field(ix,iy,iz) there is no need to update
ghost cells, whilst if you use field(ix-1,iy,iz) you do).
\\[0.5cm]
{\Large Notes\\}
The \inlinecode{field\_bc} routine just calls the
\inlinecode{do\_field\_mpi\_with\_lengths} routine which is a more general
routine that allows ghost cell information to be exchanged for fields with
an arbitrary number of cells, rather than fields which are
(-2:nx+3,-2:ny+3,-2:nz+3). This routine is used internally in the load
balancing routine when fields with both the old and new sizes must be handled
at the same time.

\pagebreak
\begin{codedef}
SUBROUTINE processor_summation_bcs
\HRule
REAL(num), DIMENSION(-2:,-2:,-2:), INTENT(INOUT) :: field
INTEGER, INTENT(IN), OPTIONAL :: flip_direction
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/boundary.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{processor\_summation\_bcs} is a routine which is used to deal with
variables, like $\vec{j}$ or number density that should be added at boundaries
to include contributions from particles on both sides of a processor boundary.
The routine is used for the current inside the particle pusher and inside most
of the routines for calculating derived variables. If you have a variable
which needs to add contributions from adjacent processors then you should
calculate the quantity on each processor, including contributions from the
particles to the ghost cells and then call this routine. When reflecting
boundary conditions are in operation, the current in the direction crossing
the boundary needs to be flipped over. This is decided upon using the
\inlinecode{flip\_direction} parameter.
\\[0.5cm]
{\Large Notes\\}
\pagebreak

These routines can be used for most MPI calls required by all but the most
extreme modifications to {\EPOCH}.

\subsection{Particle List control functions}
Collections of particles in {\EPOCH} are represented by the #particle_list#
object. These objects abstract much of the operation of the linked lists,
including adding and removing particles and sending particles to other
processors. The functions are as follows:

\pagebreak
\begin{codedef}
SUBROUTINE create_empty_partlist(partlist)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{create\_empty\_partlist} is a routine which takes a particle\_list
object and sets it up so that it points to no particles at all. It should be
used on newly allocated particle\_list objects and when a particle\_list has
served its purpose. It DOES NOT destroy the particles linked in the list at
the point that it is called. If the user wishes to delete all the particles in a
particle\_list then the routine \inlinecode{destroy\_partlist} should be used
instead.
\\[0.5cm]
{\Large Notes\\}

\pagebreak
\begin{codedef}
SUBROUTINE create_unsafe_partlist(partlist, a_particle, &
    n_elements)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
TYPE(particle), POINTER :: a_particle
INTEGER(KIND=8), INTENT(IN) :: n_elements
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{create\_unsafe\_partlist} is a routine which allows the creation of
a particle\_list which represents a subset of another particle list. This subset
is defined as starting at the particle pointed to by \inlinecode{a\_particle}
and extending for \inlinecode{n\_elements} elements. The new particle\_list is
then flagged as ``unsafe'' because if it is destroyed for any reason then it
will affect other particle lists. Many particle\_list functions can only work
on safe particle lists.
\\[0.5cm]
{\Large Notes\\}

\pagebreak
\begin{codedef}
SUBROUTINE create_unsafe_partlist_by_tail(partlist, head, &
    tail)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
TYPE(particle), POINTER :: head, tail
\end{codedef}
\vspace{1cm}
{\Large Description\\}
\inlinecode{create\_unsafe\_partlist\_by\_tail} is almost identical to
\inlinecode{create\_unsafe\_partlist}, but instead of specifying the first
particle and a number of elements, the user specifies the first and last
elements in the subset of the particle list. If the particle objects specified
for head and tail are not in the same partlist or tail actually comes before
head then the routine will fail in an undefined manner.
\\[0.5cm]
{\Large Notes\\}

\pagebreak
\begin{codedef}
SUBROUTINE create_allocated_partlist(partlist, n_elements)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
INTEGER(KIND=8), INTENT(IN) :: n_elements
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{create\_allocated\_partlist} is a helper routine to setup a new
particle\_list and create \inlinecode{n\_elements} new particle objects already
in place in the list.
\\[0.5cm]
{\Large Notes\\}
You should always use this routine when creating large numbers of new particle
objects since there is no guarantee that the internal structure of the
particle\_list objects will not change in the future. This routine will
be modified to reflect any changes in the underlying code.

\pagebreak
\begin{codedef}
SUBROUTINE create_filled_partlist(partlist, data_in, &
    n_elements)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
REAL(num), DIMENSION(:), INTENT(IN) :: data_in
INTEGER(KIND=8), INTENT(IN) :: n_elements
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{create\_filled\_partlist} is a helper routine to setup a new
particle\_list and create \inlinecode{n\_elements} new particle objects already
in place in the list. These new particle objects are then assigned properties
from the array \inlinecode{data\_in} where the particle properties are contained
in packed form. The particle data is unpacked from the array using the
\inlinecode{unpack\_particle} routine.
\\[0.5cm]
{\Large Notes\\}
You should always use this routine, if possible, when copying particles out of
packed format since there is no guarantee that the internal structure of the
particle\_list objects will not change in the future. This routine will
be modified to reflect any changes in the underlying code.

\pagebreak
\begin{codedef}
FUNCTION test_partlist(partlist)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{test\_partlist} is a routine which tests for various possible types
of error within a particle\_list object. It has a number of possible return
codes for different errors, using negative values for errors so severe that
the main tests cannot be run, or with a bitmask for errors in the main tests.
The return codes are:
\begin{itemize}
\item 0 - No error, particle\_list has passed all tests.
\item -1 - Either the head or tail of the particle\_list object is NULL. This is
  a serious error and usually means that there is a serious error inside the
  particle\_list routines.
\end{itemize}
The other error codes are returned as a bitmask and mean the following
\begin{itemize}
\item 1 - A particle\_list marked as safe has a head element which is linked to
  a preceding particle object.
\item 2 - A particle\_list marked as safe has a tail element which is linked to
  a proceding particle object.
\item 4 - The count property of a particle\_list does not correspond to the
  actual number of objects linked between the head and tail objects. This error
  code on its own usually means that the count property has been modified
  improperly.
\end{itemize}
\vspace{0.5cm}
{\Large Notes\\}
Note that this routine is only intended for debugging and is very slow. It
should never be used by the code in normal operation and all routines should be
written in such a way that it is impossible for a particle\_list object to
become corrupted.

\pagebreak
\begin{codedef}
SUBROUTINE destroy_partlist(partlist)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{destroy\_partlist} is a helper routine to delete all the particles
attached to a particle\_list and free up the memory that they use. It also
guarantees to leave the particle\_list object itself in a blank state where new
particles can be added to it. It DOES NOT delete the particle\_list object
itself, since it does not know whether or not the particle\_list is dynamically
allocated. If using dynamically allocated particle\_list objects then it is up
to the user to deallocate them AFTER the attached particles are destroyed using
\inlinecode{destroy\_partlist}.
\\[0.5cm]
{\Large Notes\\}
If a particle\_list is deleted without deleting the attached particle objects,
either using this routine or explicitly by the user, then the particles will
become orphaned and sit around using memory until the code ends. If this
happens regularly then the code will quickly crash, usually with a SIG\_SEGV
error.

\pagebreak
\begin{codedef}
SUBROUTINE copy_partlist(partlist1, partlist2)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist1, partlist2
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{copy\_partlist} is a routine which sets \inlinecode{partlist2} to
point to the same linked list of particles as \inlinecode{partlist1}. It does
not copy the particles, just sets the head and tail pointers of
\inlinecode{partlist1} to point to the same particle objects as
\inlinecode{partlist1}.
\\[0.5cm]
{\Large Notes\\}

\pagebreak
\begin{codedef}
SUBROUTINE append_partlist(head, tail)
\HRule
TYPE(particle_list), INTENT(INOUT) :: head, tail
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{append\_partlist} is a routine which takes the particles
attached to the particle\_list object \inlinecode{tail} and adds them to the
end of the linked list for particle\_list \inlinecode{head}. The particle\_list
\inlinecode{tail} is then set to be an empty particle\_list.
\\[0.5cm]
{\Large Notes\\}
This routine can only append one safe particle\_list to another safe
particle\_list.

\pagebreak
\begin{codedef}
SUBROUTINE add_particle_to_partlist(partlist, new_particle)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
TYPE(particle), POINTER :: new_particle
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{add\_particle\_to\_partlist} adds a new
particle (\inlinecode{new\_particle}) to the end of the linked list of
particles in the particle\_list object \inlinecode{partlist}. It deals with
cases of empty particle\_list objects automatically.
\\[0.5cm]
{\Large Notes\\}
If you want to add a new particle to the end of a particle list you should
always use this routine.

\pagebreak
\begin{codedef}
SUBROUTINE remove_particle_from_partlist(partlist, &
    a_particle)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
TYPE(particle), POINTER :: a_particle
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{remove\_particle\_from\_partlist} removes the particle object
specified by \inlinecode{a\_particle} from the particle\_list object given by
\inlinecode{partlist}. Be very careful that \inlinecode{a\_particle} is indeed
in the linked list pointed to by \inlinecode{partlist}, otherwise it is possible
for the particle\_list object which really does contain \inlinecode{a\_particle}
to be left with an invalid pointer as its head or tail element if
\inlinecode{a\_particle} is either the head or tail element.
\\[0.5cm]
{\Large Notes\\}
Although this routine does work with unsafe particle\_list objects, you should
be very careful using it in this case as it can break the head or tail element
of the primary particle\_list which the unsafe particle\_list is a subset of.
As a general rule, you should only use this routine to remove particles from a
simple particle\_list which is a singly referenced primary, safe particle\_list.

\pagebreak
\begin{codedef}
SUBROUTINE setup_partlists()
\HRule

\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{setup\_partlists} is a routine which is called once when {\EPOCH}
first starts. It sets the variable \inlinecode{nvars} which is the number of
REAL(num) values required to contain all the information about a
single particle object needed when a particle is transferred to another
processor. How the information is packed and unpacked from the particle object
into an array of REAL(num) values is controlled in the functions
\inlinecode{pack\_particle} and \inlinecode{unpack\_particle}.
\\[0.5cm]
{\Large Notes\\}
If the particle type gains additional properties as the result of preprocessor
directives then there should be a line which increments \inlinecode{nvars} by
the correct number when that preprocessor directive is active. For example:
\begin{boxverbatim}
#ifdef PER_PARTICLE_CHARGE_MASS
  nvar = nvar+2
#endif
\end{boxverbatim}

\pagebreak
\begin{codedef}
SUBROUTINE   pack_particle(array, a_particle)
SUBROUTINE unpack_particle(array, a_particle)
\HRule
REAL(num), DIMENSION(:), INTENT(INOUT) :: array
TYPE(particle), POINTER :: a_particle
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{pack\_particle} and \inlinecode{unpack\_particle} are subroutines
which are used to copy all the information about a particle
necessary for the particle to be transferred to another processor into a
temporary array before sending to another processor. If a new particle property
has been added to the particle then these routines must be modified to allow
the copying of the new data into the array. The parameter \inlinecode{array} is
a REAL(num) array of length \inlinecode{nvars} and is the array into which the
data either must be packed or from which it must be
unpacked. \inlinecode{a\_particle} is the particle object which must either have
its data copied into the array or be
populated with data from the array. No restriction is placed on how the
data should be packed into the data array, but obviously
\inlinecode{pack\_particle} and \inlinecode{unpack\_particle} must be inverse
operations so that particles packed by one processor can be unpacked correctly
by another processor.
\\[0.5cm]
{\Large Notes\\}
Since it is very unlikely that {\EPOCH} will be run on anything other than a
homogeneous cluster, it is acceptable to use the Fortran \inlinecode{TRANSFER}
function to pack incompatible data types into the \inlinecode{array} array. Just
make sure that \inlinecode{nvars} is defined in \inlinecode{setup\_partlists}
to be long enough to contain all the information. More documentation on the
\inlinecode{TRANSFER} function (which is rarely used and dangerous!) can be
found at \url{http://www.macresearch.org/%
advanced_fortran_90_callbacks_with_the_transfer_function}.

\pagebreak
\begin{codedef}
SUBROUTINE display_particle(a_particle)
\HRule
TYPE(particle), POINTER :: a_particle
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
%\setlength{\emergencystretch}{3em}
Displays the key information about a particle given by the parameter
\inlinecode{a\_particle}. Used by \inlinecode{compare\_particles}.
\\[0.5cm]

\pagebreak
\begin{codedef}
FUNCTION compare_particles(particle1, particle2)
\HRule
TYPE(particle), POINTER :: particle1, particle2
LOGICAL :: compare_particles
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
Compares all the properties of two particle objects and displays the
information if they don't match. Used internally by
\inlinecode{test\_packed\_particles}. If the particle object is extended then
this routine should also be modified to test for equivalence of the new
properties.
\\[0.5cm]

\pagebreak
\begin{codedef}
FUNCTION test_packed_particles(partlist, array, &
    npart_in_data)
\HRule
TYPE(particle_list), INTENT(IN) :: partlist
REAL(num), DIMENSION(:), INTENT(IN) :: array
INTEGER(KIND=8), INTENT(IN) :: npart_in_data
LOGICAL :: test_packed_particles
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{test\_packed\_particles} is a routine which checks that a packed
array of particles can be successfully unpacked back into particle objects. The
parameters are:
\begin{itemize}
\item \inlinecode{partlist} - The particle\_list corresponding to the original
  unpacked particles.
\item \inlinecode{array} - The REAL(num) array containing the packed data.
\item \inlinecode{npart\_in\_data} - The number of particles which were packed
  into \inlinecode{array}.
\end{itemize}
The routine tests that the number of particles in the particle\_list match the
number believed to be in the data array, that the length of the data array is
correct and then unpacks each particle in turn from the data array and uses the
\inlinecode{compare\_particles} function to compare the particles with the
original versions in the particle\_list. If any particles fail the comparison
then an error is output to stdout and the function returns a value of
#.FALSE.#. The error message includes the processor
rank on which the problem occurs but the routine does not specifically include
any MPI commands, so it is possible to call the routine on a subset of
processors.
\\[0.5cm]

\pagebreak
\begin{codedef}
SUBROUTINE partlist_send(partlist, dest)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
INTEGER, INTENT(IN) :: dest
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{partlist\_send} is a routine for sending all the particles in the
particle\_list object \inlinecode{partlist} to another processor. The
destination processor is identified by its rank which is given by the
\inlinecode{dest} parameter. The routine does not destroy the particle\_list
object which is given to it.
\\[0.5cm]
{\Large Notes\\}
\inlinecode{partlist\_send} uses MPI blocking sends, so unless a matching
\inlinecode{partlist\_recv} has been posted on \inlinecode{dest} then the
routine will deadlock. It would be fairly simple to write a non-blocking
version of \inlinecode{partlist\_send}, but at present no need for such a
routine has been found.

\pagebreak
\begin{codedef}
SUBROUTINE partlist_recv(partlist, src)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist
INTEGER, INTENT(IN) :: src
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{partlist\_recv} is a routine for receiving particles sent by a call
to \inlinecode{partlist\_send} and loading them into the particle\_list object
\inlinecode{partlist}. The source processor is identified by its rank which is
given by the \inlinecode{src} parameter. The routine destroys the particle\_list
which it is given and indeed will leave orphaned particles if it is not given
an empty particle\_list to receive the data.
\\[0.5cm]
{\Large Notes\\}
\begin{itemize}
\item \inlinecode{partlist\_recv} uses MPI blocking receives, so unless a
  matching \inlinecode{partlist\_send} has been posted on \inlinecode{src} then
  the routine will deadlock. It would be fairly simple to write a non-blocking
  version of \inlinecode{partlist\_recv}, but at present no need for such a
  routine has been found.

\item Although it is not possible to directly use \inlinecode{partlist\_recv}
  to add new particles onto an existing particle\_list, it is only two lines to
  do this. First call \inlinecode{partlist\_recv} with a temporary
  particle\_list to receive the data and then use \inlinecode{append\_partlist}
  to attach the particles to the end of the already populated list.
\end{itemize}

\pagebreak
\begin{codedef}
SUBROUTINE partlist_send_recv(partlist_send, &
    partlist_recv, dest, src)
\HRule
TYPE(particle_list), INTENT(INOUT) :: partlist_send
TYPE(particle_list), INTENT(INOUT) :: partlist_recv
INTEGER, INTENT(IN) :: dest, src
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/partlist.F90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{partlist\_sendrecv} is a routine equivalent to
\inlinecode{MPI\_SENDRECV} in that it allows overlapping sends and receives to
be written in a single line rather than the end user having to split processors
into red/black ordered pairs for communication. It sends the particle data in
\inlinecode{partlist\_send} to the processor with rank \inlinecode{dest} and
receives particle data sent by processor \inlinecode{src} and stores it in the
particle\_list \inlinecode{partlist\_recv}. The routine is destructive to both
sending and receiving particle\_lists, and can lead to orphaned particles if a
filled particle\_list is passed as \inlinecode{partlist\_recv}. this is the
routine which is used in the particle boundary conditions.
\\[0.5cm]
{\Large Notes\\}
\begin{itemize}
\item \inlinecode{partlist\_sendrecv} uses MPI blocking sendrecv commands, so
  should be used in matching pairs or the routine will deadlock.

\item Although it is not possible to directly use
  \inlinecode{partlist\_sendrecv} to add new particles onto an existing
  particle\_list, it is only two lines to do this. First call
  \inlinecode{partlist\_sendrecv} with a temporary particle\_list to receive the
  data and then use \inlinecode{append\_partlist} to attach the particles to
  the end of the already populated list.
\end{itemize}
\pagebreak


\subsection{MPI in {\EPOCH}}
{\EPOCH} is a massively parallel code written using standard MPI1.2, and likely
to be upgraded to MPI2 in the near future. Due to the massively parallel nature
of {\EPOCH}, there are MPI commands scattered throughout many parts of the code,
although the MPI has been hidden as far as possible from the end user. The main
use of MPI occurs during I/O, in the boundary conditions and during load
balancing. The MPI setup routines are all in
\inlinecode{src/housekeeping/mpi\_routines.F90}, and the routines which are
used to create the MPI types used by MPI-IO are in
\inlinecode{src/housekeeping/mpi\_subtype\_control.f90}.

\subsubsection{General MPI in {\EPOCH}}
{\EPOCH} uses Cartesian domain decomposition for parallelism and creates an MPI
Cartesian topology using \inlinecode{MPI\_CART\_CREATE}. The use of MPI in
{\EPOCH} is deliberately kept as simple as possible, but there are some points
which must be made and some variables which must be explained.
\begin{itemize}
\item MPI decomposition is reversed compared to array ordering. Due to the
  layout of arrays in Fortran, you get slightly faster performance if you split
  arrays so that the first index remains as long as possible. Since {\EPOCH}
  uses \inlinecode{MPI\_DIMS\_CREATE} to do array subdivision, this means that
  the MPI coordinate system is ordered backwards compared to the main arrays.
  This means that the \inlinecode{coordinates} array which holds the
  coordinates of the current processor in the Cartesian topology is ordered
  as \{coord\_z, coord\_y, coord\_x\}.
\item To make this easier, there are some helper variables. The simplest of
  these just gives the processors attached to each face of the domain on the
  current processor. These variables are named \inlinecode{x\_min, x\_max,
  y\_min, y\_max, z\_min} and \inlinecode{z\_max}.
\item Since it is possible for particles to cross boundaries diagonally there
  is another variable \inlinecode{neighbour} which identifies every possible
  neighbouring processor including those meeting at single edges and at
  corners. \inlinecode{neighbour} is an array which runs (-1:1,-1:1,-1:1) and,
  perhaps inconsistently, is ordered in normal order rather than reversed
  order. This means that \inlinecode{x\_min == neighbour(-1,0,0)} and
  \inlinecode{z\_max == neighbour(0,0,1)}.
\item The variable \inlinecode{comm} is the handle for the Cartesian
  communicator returned from \linebreak MPI\_CART\_CREATE.
\item The variable \inlinecode{errcode} is the standard error variable for all
  MPI communications. However, {\EPOCH} uses the standard
  MPI\_ERRORS\_ARE\_FATAL error handler so this variable is never tested.
\item {\EPOCH} uses a single variable, \inlinecode{status}, to hold all MPI
  status calls. Since there is no non-blocking communication this variable
  is never checked.
\item The rank of the current processor is stored in the variable
  \inlinecode{rank}.
\item The number of processors is stored in \inlinecode{nproc}.
\item The number of processors assigned to any given direction of the Cartesian
  topology is given by \inlinecode{nproc\{x,y,z\}}.
\end{itemize}

There are some other variables which are not technically part of the MPI
implementation, but which only exist because the code is running in
parallel. These are
\begin{itemize}
\item #REAL(num) :: {x,y,z}_min_local# - The location of the start of the
  domain on the local processor in real units.
\item #REAL(num) :: {x,y,z}_max_local# - The location of the end of the
  domain on the local processor in real units.
\item #INTEGER, DIMENSION(1:nproc{x,y,z}) :: cell_{x,y,z}_min# - The cell
  number for the start of the local part of the global array in each direction.
\item #INTEGER, DIMENSION(1:nproc{x,y,z}) :: cell_{x,y,z}_max#
  - The cell number for the end of the local part of the global array in each
  direction.
\end{itemize}

\subsubsection{\inlinecode{mpi\_routines.F90}}
\inlinecode{mpi\_routines.F90} is the file which contains all the MPI setup
code. It contains the following routines:
\begin{itemize}
\item mpi\_minimal\_init - Contains code to start MPI enough to
  allow the input deck reader to work. The default {\EPOCH} code setup means
  that it needs to initialise MPI, obtain the rank and the number of processors.
\item setup\_communicator - Routine which creates the Cartesian communicator
  used by the code after the input deck has been parsed. It also populates
  \inlinecode{x\_min, x\_max} etc. It is in its own subroutine so that it can be
  recalled after the start of the window move when the code is using a moving
  window. This is needed since it is valid to have a non-periodic boundary
  before the window starts to move and a periodic boundary afterwards.
\item mpi\_initialise - This routine calls \inlinecode{setup\_communicator} and
  then allocates all the arrays to do with fields, etc. It also sets up the
  particle list objects for each species. If the code is running with only
  manual initial conditions then this routine loads the requested number of
  particles on each processor. Otherwise either the restart or the autoloader
  code load the particles.
\item mpi\_close - This routine performs all the needed cleanup before the
  final call to \inlinecode{MPI\_FINALIZE}.
\end{itemize}

\subsubsection{\inlinecode{mpi\_subtype\_control.f90}}
This file contains all the routines which are used to create the MPI types
which are used in the SDF I/O system. Most of the routines in this section are
used to create the types used for writing the default variables and,
when modifying the code, it is possible to output anything which has the same
shape and size on disk as the default variables without ever having to use the
routines in this file. However, if you are creating more general modifications
which can include variables of different sizes with different layouts across
processors then you may wish to use these routines to create new MPI types which
match your data layout. Any valid MPI type describing the data layout will work
with the SDF library, so there is no absolute need to use these routines. Only
the general purpose subroutines are described here, since most of the other
routines are fairly clear and use these routines internally.

\pagebreak
\begin{codedef}
FUNCTION create_particle_subtype(npart_local)
\HRule
INTEGER(KIND=8), INTENT(IN) :: npart_local
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/mpi\_subtype\_control.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{create\_particle\_subtype} is a routine which creates an MPI type
representing particles which are spread across different processors with
\inlinecode{npart\_local} particles on each
processor. \inlinecode{npart\_local} does not have to be the same number on all
processors.
\\[0.5cm]
{\Large Notes\\}
Currently this is only used for reading particle data from restart snapshots.
It is likely to go away in the near future.

\pagebreak
\begin{codedef}
SUBROUTINE create_ordered_particle_offsets(n_dump_species,&
    npart_local)
\HRule
INTEGER, INTENT(IN) :: n_dump_species
INTEGER(KIND=8), DIMENSION(n_dump_species), &
    INTENT(IN) :: npart_local
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/mpi\_subtype\_control.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{create\_ordered\_particle\_offsets} is a routine which creates an
array of offsets representing particles from \inlinecode{n\_dump\_species}
which are spread across different processors with
\inlinecode{npart\_local(ispecies)} particles of each species on each
processor. \inlinecode{npart\_local} does not
have to be the same number on all processors and does not have to be the same
number for each species.
\\[0.5cm]
{\Large Notes\\}

\pagebreak
\begin{codedef}
FUNCTION create_field_subtype(n\{x,y,z\}_local, &
    cell_start_\{x,y,z\}_local)
\HRule
INTEGER, INTENT(IN) :: nx_local, ny_local, nz_local
INTEGER, INTENT(IN) :: cell_start_x_local
INTEGER, INTENT(IN) :: cell_start_y_local
INTEGER, INTENT(IN) :: cell_start_z_local
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/housekeeping/mpi\_subtype\_control.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{create\_field\_subtype} is a routine which creates an MPI type
representing a field that is defined across some or all of the processors. The
\inlinecode{n\{x,y,z\}\_local} parameters are the number of points in the x,
y, z directions (if they exist in the version of the code that you are working
on) that are on the local processor. The
\inlinecode{cell\_start\_\{x,y,z\}\_local} parameters are the offset of the
top, left, back corner of the local subarray in the global array that would
exist if the code was running on one processor. This is an {\it offset}, not a
position and so it begins at \{0,0,0\} NOT \{1,1,1\}.
\\[0.5cm]
{\Large Notes\\}
In EPOCH3D there is also a routine called
\inlinecode{create\_field\_subtype\_2d} which is exactly equivalent to
\inlinecode{create\_field\_subtype} in EPOCH2D and is used for writing the 2D
distribution functions. At present, there are not equivalent 1D functions
except in EPOCH1D, but these could easily by added if required.

\pagebreak

\subsection{The load balancer}
One of the major limiting factors in the scalability of PIC codes is load
balancing. Due to the synchronisation of the currents required for the update
of the EM fields the entire code runs at the speed of the slowest
process. Since most of the time in the main {\EPOCH} cycle is taken by the
particle pusher, this equates to the process with the highest number of
particles being the slowest. Since the location of particles is dependent upon
the solution of the problem under consideration, in general the code will not
have exactly the same number of particles on each processor. The load balancer
is used to move the inter-processor boundaries so that the number of particles
is as close to the same on each processor as possible. The load balancer is
invoked at the start of the code and when the ratio of the least loaded
processor to the most loaded processor falls below a user specified critical
point.\\

{\EPOCH}'s load balancer works by rearranging the processor boundaries in 1D
sweeps in each direction, rather than attempting to perform multidimensional
optimisation. Also, at present the MPI in {\EPOCH} requires each processor to be
simply connected at every point, so it must have one processor to the left, one
to the front etc. which introduces a further restriction on the load
balancer. Otherwise, the load balancer is fully general. The load balancing
sweep is illustrated in Figure~\ref{sweep}.
The load balancer is
implemented in the file \inlinecode{src/housekeeping/balance.F90} and is called
by the routine \inlinecode{balance\_workload(override)}. The parameter
\inlinecode{override} is used to force the code to perform a load balancing
sweep even when it would normally determine that the imbalance is not large
enough to force a load balancing sweep. Although the load balancer is hard
coded to load balance in all available directions, the code is written in such
a way that it is possible to modify the code to load balance in only one
direction, or to automatically determine which single direction gives the best
performance.

\captionedimage{./images/sweep}{sweep}{Illustration of the load balancing sweep}

The details of the load balancer are fairly intricate, and if major
modification to the load balancer is required, it is recommended that the
original authors be contacted for detailed advice on how to proceed. However,
the general layout of the routine is as follows.

\begin{itemize}
\item Use MPI\_ALLREDUCE to determine the global minimum and maximum number of
  particles. If the ratio of the minimum to the maximum is above the load
  balance threshold then just return from the subroutine.
\item The code uses the routines \inlinecode{get\_load\_in\_\{x,y,z\}} to
  determine the work load along each direction of the domain.
  The \inlinecode{get\_load\_in\_x} routine uses the x-coordinate of every
  particle to create a 1D particle density in the x-direction. This is then
  combined with the total number of grid cells in the y-z plane to give
  a 1D array of the work load in the x-direction. Similarly for
  \inlinecode{get\_load\_in\_\{y,z\}}.
\item Next the load array is passed to the \inlinecode{calculate\_breaks}
  routine which fills the arrays \inlinecode{starts\_\{x,y,z\}}
  and \inlinecode{ends\_\{x,y,z\}}. These arrays contain the starting and
  ending cell numbers of the hypothetical global array in each direction for
  each processor.
\item The routine \inlinecode{redistribute\_fields} is then called to move the
  information about field variables which cannot be recalculated. If new field
  variables are created that cannot be recalculated after the load balancing is
  completed then \inlinecode{redistribute\_fields} has to be modified for these
  new variables.
\item The next section of the routine deals with those variables which can be
  recalculated after the load balance sweep is complete, such as the coordinate
  axes and the arrays which hold the particle kinetic energy.
\item The penultimate section of the routine then changes the variables which
  tell the code where the edges of its domain lie in real space to reflect the
  changed shape of the domains.
\item The final part is the call to \inlinecode{distribute\_particles} which
  moves the particles to the new processor. Once this is
  finished, the code should have as near as possible the same number of
  particles on each processor.
\end{itemize}

Most of the load balancer is purely mechanical and should only be changed if
the way in which the code is to perform load balancing is fundamentally
altered. The redistribution of particles that takes place in
\inlinecode{distribute\_particles} uses the standard particle\_list objects, so
that if the necessary changes have been made to the routines in
\inlinecode{src/housekeeping/partlist.F90} to allow correct
boundary swaps of particles then the load balancer should work with no further
modification. The only part of the load balancer which should need changing is
\inlinecode{redistribute\_fields} which requires explicit modification if new
field variables are required. For fields which are the same shape as the main
array there is significant assistance provided within the code to make the
re-balancing simpler. There are also routines which can help with re-balancing
variables which are the size of only one edge or face of the domain. Variables
which are of a completely different size but still need to be rebalanced when
coordinate axes move have to have full load balancing routines implemented by
the developer. This is beyond the scope of this manual and any developer who
needs assistance with development such a modification should contact the
original author of the code. The field balancer is fairly simple and mostly
calls one of three routines: \inlinecode{redistribute\_field} and either
\inlinecode{redistribute\_field\_2d} or \inlinecode{redistribute\_field\_1d}
depending on the dimensionality of your code. To redistribute full field
variables the routine to use is \inlinecode{redistribute\_field}, and an
example of using the code looks like:
\begin{boxverbatim}
  temp = 0.0_num
  CALL redistribute_field(new_domain, bz, temp)
  DEALLOCATE(bz)
  ALLOCATE(bz(-2:nx_new+3,-2:ny_new+3))
  bz = temp
\end{boxverbatim}

In this calling sequence the
\inlinecode{redistribute\_field} subroutine is used to redistribute the field
\inlinecode{bz}, and the newly redistributed field is copied
into \inlinecode{temp}; an array which is already allocated to the
correct size. The \inlinecode{new\_domain} parameter is an array indicating the
location of the start and end points of the new domain for the current processor
in gridpoints offset from the start of the global array. It is passed into the
\inlinecode{redistribute\_fields} subroutine as a parameter from the
\inlinecode{balance\_workload} subroutine and should not be changed. The
\inlinecode{temp} variable is needed since Fortran standards before Fortran2000
do not allow the deallocation and reallocation of parameters passed to a
subroutine. There is a more elegant solution, where \inlinecode{temp} is
hidden inside the \inlinecode{redistribute\_field} subroutine. However, support
for this in current Fortran2000 implementations is unreliable.\\

The routine for re-balancing variables which lie along an edge of the domain are
very similar and are demonstrated in the \inlinecode{redistribute\_fields}
subroutine for lasers attached to different boundaries. It is
recommended that a developer examine this code when developing new routines.

\subsection{{\EPOCH} I/O}
\label{sec:io}

{\EPOCH} uses a file format called {\it SDF} which was custom developed for use
with codes developed by CFSA at the University of Warwick. The internal
structure of this format is documented in \sect{sdf},
but it is not necessary to have a full understanding of the file format
to add the output of new variables to {\EPOCH}. To add a new variable to
{\EPOCH}'s output, you simply have to use the supplied subroutines of the SDF
library which is part of {\EPOCH}. The file output from {\EPOCH} takes place in
{\it diagnostics.F90}, so to add new variables to the output you must
add additional code there. Looking through the listings, you will see two lines:
\begin{boxverbatim}
  CALL sdf_open(sdf_handle, filename, rank, comm, c_sdf_write)

  CALL sdf_close(sdf_handle)
\end{boxverbatim}

These, as may be expected, are the commands which open and close the SDF file.
It is perfectly possible to create new SDF files containing only your own data.
There are various commands in-between which actually write the data into the
file. These commands start with \inlinecode{sdf\_} to ensure that the don't
conflict with any other subroutine names in the code.
Some more complex areas of I/O, such as the particle probes
and the distribution function routines call other subroutines in their
respective source files, but these too make use of the SDF routines to actually
write data. A user should never try to write data directly to the output file,
since this will cause problems with internal parts of the SDF format and
generate a nonsensical file.

\subsubsection{dumpmask}
Looking through {\it diagnostics.F90} there are many lines with commands which
begin \inlinecode{sdf\_}, but are all prepended with a command which looks
like:
\begin{boxverbatim}
  IF (IAND(dumpmask(c_dump_id), code) .NE. 0) THEN
\end{boxverbatim}
This is the method by which {\EPOCH} allows the end user to specify whether a
variable should be dumped, and whether it should only be dumped at
full/partial/restart dumps. #dumpmask# is an integer array, the length of
which is defined by the variable \inlinecode{num\_vars\_to\_dump} in
\inlinecode{shared\_data.F90} and contains the bitmask representing all the
types of output which should be written for the associated variable. The
possible values in the bitmask are:

\begin{itemize}
\item \inlinecode{c\_io\_never} - Never dump this variable.
\item \inlinecode{c\_io\_always} - Dump this variable at every output dump.
\item \inlinecode{c\_io\_full} - Dump this variable at full dumps.
\item \inlinecode{c\_io\_restartable} - Dump this variable for restart dumps.
\item \inlinecode{c\_io\_species} - If meaningful for this variable, write
  information for each species rather than integrated over all species.
\item \inlinecode{c\_io\_no\_sum} - If meaningful for this variable, do not
  write information integrated over all species.
\item \inlinecode{c\_io\_averaged} - If meaningful for this variable, write
  this variable averaged over time.
\item \inlinecode{c\_io\_snapshot} - If meaningful for this variable, write
  the non-averaged value of the variable.
\end{itemize}

The \inlinecode{c\_dump\_id} entry is a constant defined in
\inlinecode{shared\_data.F90} which identifies the variable's index within
this dumpmask.

When adding a new variable to be written to disk, the value of
\inlinecode{num\_vars\_to\_dump} should be increased to match the number of new
written variables. Next, open the file \inlinecode{src/deck/deck\_io\_block.F90}
and find the line:
\begin{boxverbatim}
  CHARACTER(LEN=entry_length), DIMENSION(io_block_elements) :: &
      io_block_name = ...
\end{boxverbatim}

Simply add new strings for your new variables to the end of the
definitions along with its \inlineemph{c\_dump\_id} value. These new variable
names should then be placed in your input decks
in the same place as the existing I/O information and take the same parameters.

\subsection{Adding a derived variable}

As already stated, derived variables are variables which are defined on the
Cartesian spatial grid but are not directly updated by the solver. They are
calculated when needed for output or for use in other physics packages. The
final form of a derived variable is an array on each processor with the same
size as the field arrays.

\subsubsection{Writing your variable to disk}
A derived variable is defined on the same grid as the main simulation variables
and must be written to disk in such a way as to stitch the parts of the grid
from each processor together. This is achieved using the routine:
\begin{boxverbatim}
CALL sdf_write_plain_variable(sdf_handle, id, name, units, dims, stagger, &
    grid_id, variable, subtype_field, subarray_field)
\end{boxverbatim}
The parameters have the following types and meanings:

\begin{itemize}
\item block\_id - The id name of the variable. This character string is a
  unique identifier for
  the variable in the file enabling a program to retrieve it later. Once
  defined it should not change so that newer versions of {\EPOCH} can still
  identify variables generated by older versions.
\item name - The display name of the variable. This character string is
  the name that is used
  by external programs to display an identifying name for the variable. If it
  contains '/' characters then these are used by VisIt to group the variables.
\item units - The units of the variable. This character string is used when
  displaying the data units. For most variables in {\EPOCH} these
  are SI units.
\item dims - An nD integer array containing the GLOBAL length of the variable
  across all processors. In {\EPOCH} a variable actually called ``dims'' exists
  for variables which are the same size as the default field variables.
\item stagger - An integer constant containing the stagger of a variable from
  the cell centre of a cell. This property lets external programs know the
  position of a variable on the grid.
\item grid\_id - The id name of the grid to which the variable is attached. In
  {\EPOCH}, the main grid is just called ``grid''. Note that this property is
  case sensitive.
\item variable - The actual variable to be written to disk.
\item subtype\_field - This is an MPI type representing the layout of the data
  across the processors. For a standard field variable, there is an
  automatically created type called ``subtype\_field'' which should be used
  here.
\item subarray\_field - This is an MPI type representing the section of
  the ``variable'' parameter to be written. For a standard field variable,
  there is an automatically created type called ``subarray\_field'' which
  should be used here.
\end{itemize}

It's probably easiest to read the diagnostics.F90 file and see how the code
implements the output of simple variables like ex or ey for an example of how
this works. Once the appropriate sdf\_write call has
been added to the code, there is no further work
to be done. The IDL, MatLab and VisIt routines will all read the existence of
the variable from the metadata in the SDF file, and it will now be available to
view in all SDF reading packages.

There is a working variable called \inlinecode{array} which is large enough to
store a derived variable. It is
therefore recommended that to calculate derived variables a new subroutine
should be created which populates \inlinecode{array} with the required variable
and then writes it to disk. An example would look like:
\begin{boxverbatim}
IF (IAND(dumpmask(c_dump_myvar), code)) THEN
  CALL calc_my_variable(array)

  CALL sdf_write_plain_variable(sdf_handle, 'my_var', 'Mine/variable', 'unit',
      dims, c_stagger_cell_centre, 'grid', array, subtype_field, subarray_field)
ENDIF
\end{boxverbatim}
where \inlinecode{calc\_my\_variable} is a function which calculates the
variable which you wish to write. The form of this function depends on the type
of variable to be calculated and is given in the next section.

\subsection{Adding a particle variable}
The next simplest type of output to add is a new property for all particles. To
add new particle variables to the output dump, two things are needed: a call
to the SDF command to write the data and an iterator function to iterate
through all the particles. The iterators are stored in the file
{\it iterators.F90}. An example iterator is:
\begin{boxverbatim}
  ! iterator for particle momenta
  FUNCTION iterate_px(array, n_points, start)

    REAL(num) :: iterate_px
    REAL(num), DIMENSION(:), INTENT(OUT) :: array
    INTEGER, INTENT(INOUT) :: n_points
    LOGICAL, INTENT(IN) :: start
    TYPE(particle), POINTER, SAVE :: cur
    TYPE(particle_list), POINTER, SAVE :: current_list
    INTEGER :: part_count

    IF (start)  THEN
      CALL start_particle_list(current_species, current_list, cur)
    ENDIF

    part_count = 0
    DO WHILE (ASSOCIATED(current_list) .AND. (part_count .LT. n_points))
      DO WHILE (ASSOCIATED(cur) .AND. (part_count .LT. n_points))
        part_count = part_count + 1
        array(part_count) = cur%part_p(1)
        cur=>cur%next
      ENDDO
      ! If the current partlist is exhausted, switch to the next one
      IF (.NOT. ASSOCIATED(cur)) CALL advance_particle_list(current_list, cur)
    ENDDO
    n_points = part_count

    iterate_px = 0

  END FUNCTION iterate_px
\end{boxverbatim}

This is a fairly complicated routine which includes code for dealing with the
possibility of particle species not being dumped, and other complicated
book keeping. Luckily, there is only one line in the routine which needs to
change to output a new variable. The is the line:
\begin{boxverbatim}
        array(part_count) = cur%part_p(1)
\end{boxverbatim}

To write a new iterator, you just have to copy the skeleton of an existing
iterator and change this line to copy your particle property into the ``array''
array. The details of the particle structure's contents is explained in
\sect{partrep}. Once your new iterator has been written and added into
the {\it iterators.F90} file, it's time to add the SDF routine to actually
write the data. The routine is:
\begin{boxverbatim}
    CALL write_particle_variable(c_dump_id, code, name, iterator)
\end{boxverbatim}

The parameters this time are
\begin{itemize}
\item c\_dump\_id - The index into the dumpmask for this variable.
\item code - The dump code for the current output dump.
\item name - The display name to use for this variable.
\item iterator - The name of the iterator function that you created in
  the previous step. Note that this is not a string but simply the name of the
  function.
\end{itemize}

Once again, looking at how this is implemented for one of the existing
variables (e.g. px) is probably the most enlightening way to see how it
works. As for the fluid variables, the new variable will appear in IDL, MatLab
and VisIt.

At this point it is possible to write any property which is similar to the
default field variables or the default particle properties. It becomes slightly
more challenging if you want to write other types of variable into an output
file.

\pagebreak

\subsection{Precompiler directives}
{\EPOCH} uses precompiler directives to switch certain features of the code on
or off. The precompiler directives all begin with a ``\#'' character and look
like:
\begin{boxverbatim}
#ifdef MY_PRECOMPILER_DIRECTIVE
  some_fortran_of_some_kind
#else
  some_other_fortran
#endif
\end{boxverbatim}
They behave in a very simple manner. The precompiler runs BEFORE the
Fortran compiler and, until it reaches a precompiler directive, it just creates
a temporary file which is an exact copy of the source file. When it reaches a
precompiler directive of this kind it treats the \#ifdef commands as
if/then/else statements. If
\inlinecode{MY\_PRECOMPILER\_DIRECTIVE} was defined in the makefile then
\inlinecode{some\_fortran\_of\_some\_kind} is pushed out to the temporary
file. Otherwise \inlinecode{some\_other\_fortran} is written instead.
The precompiler directives themselves are never output to the temporary
file. Once then preprocessor has finished, it passes this temporary file
to the Fortran compiler which can then compile it just like any other
standard Fortran file.

\subsubsection{When to use precompiler directives}
\begin{itemize}
\item When adding properties to the \inlinecode{particle} structure.
\item When adding time consuming calculations to the particle pusher.
\end{itemize}
Precompiler directives should be avoided when there is no significant
performance gain or memory reduction to be made. Wherever possible, optional
features should be controlled by parameters in the input deck.

\subsubsection{The directive printing routine on code startup}
When {\EPOCH} starts it prints the precompiler directives that it was built with
and what they mean. This isn't required, but has proved very useful and is
implemented in a very simple way. Just open the file
\inlinecode{src/housekeeping/welcome.F90} and find the subroutine
\inlinecode{compiler\_directives}. There are a large block of precompiler
directives which read:
\begin{boxverbatim}
#ifdef TRACER_PARTICLES
    defines = IOR(defines, c_def_tracer_particles)
    WRITE(*, *) "Tracer particle support -DTRACER_PARTICLES"
#endif
\end{boxverbatim}

Simply add a new element to the end of the list.
\begin{boxverbatim}
#ifdef MY_PRECOMPILER_DIRECTIVE
    defines = IOR(defines, c_def_my_precompiler_directive)
    WRITE(*,*) "My new physics -DMY_PRECOMPILER_DIRECTIVE"
#endif
\end{boxverbatim}

You will also need to add #c_def_my_precompiler_directive# to the list of
constants in #src/shared_data.F90#.

\subsubsection{Precompiler directives and the input deck}
In theory, it is possible for someone to request a feature of the code in the
input deck which this version hasn't been compiled with. In this
case, there is a special error code \inlinecode{c\_err\_pp\_options\_wrong}
which causes the input deck parser to give a meaningful error. You should also
set the string \inlinecode{extended\_error\_string} to be the define
command for the missing preprocessor directive i.e
\inlinecode{extended\_error\_string = "-DMY\_PRECOMPILER\_DIRECTIVE"}

\section{{\EPOCH} front end programming}

\subsection{Strings in {\EPOCH}}
Fortran is not a language famous for its string handling capabilities, but due
to the presence of the input deck {\EPOCH} has fairly extensive string handling
routines. Strings used are all of the standard Fortran \inlinecode{CHARACTER}
type and are defined as:
\begin{boxverbatim}
CHARACTER(LEN=entry_length) :: string
\end{boxverbatim}
\inlinecode{entry\_length} is a global constant defined in
\inlinecode{src/shared\_data.F90} which can be increased to allow {\EPOCH} to
handle longer strings. There may be reasons to increase this length if you wish
to use long complex expressions in the input deck. Note that many Fortran
compilers do not allow strings to exceed 512 characters in length.

\subsubsection{{\EPOCH} string handling routines}

Listed here are the string handling routines (other than those in the core
maths parser routine which are documented elsewhere) which are currently used
in {\EPOCH}.

\pagebreak
\begin{codedef}
FUNCTION str_cmp(str_in, str_test)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in, str_test
LOGICAL :: str_cmp
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{str\_cmp} is the routine which does all the string comparisons in
{\EPOCH}. It deals with leading and trailing whitespace automatically and
tests for length differences. It does not test for strings being
valid substrings of each other, only for full equality.
\\[0.5cm]
{\Large Notes\\}
A developer should always use \inlinecode{str\_cmp} rather than doing their own
string testing to ensure consistent behaviour across the entire {\EPOCH} code
base.

\pagebreak
\begin{codedef}
FUNCTION as_real_simple(str_in, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
REAL(num) :: as_real_simple
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{as\_real\_simple} is a routine to convert a string into a real
number without invoking the maths parser. It can cope with standard form as
well as simple decimal reals. It is significantly faster than the maths parser,
but should only be used when the user explicitly {\it shouldn't} be able to use
a mathematical expression. If the string cannot be parsed then the routine sets
the bitmask \inlinecode{c\_err\_bad\_value} on the parameter \inlinecode{err}
\\[0.5cm]
{\Large Notes\\}
If you have a string which has to be converted into a real quickly then this is
the routine to use. You probably shouldn't use it when parsing a string from
the input deck, since there is no reason to restrict the user from specifying a
mathematical expression. The routine is used inside the maths parser to parse
simple numbers.

\pagebreak
\begin{codedef}
FUNCTION as_integer_simple(str_in, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
INTEGER :: as_integer_simple
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{as\_integer\_simple} is a routine to convert a string into an
integer without invoking the maths parser. It can cope with standard form as
well as simple decimal integers. It is significantly faster than the maths
parser, but should only be used when the user explicitly {\it shouldn't} be
able to use a mathematical expression. If the string cannot be parsed then the
routine sets the bitmask \inlinecode{c\_err\_bad\_value} on the parameter
\inlinecode{err}
\\[0.5cm]
{\Large Notes\\}
This routine is used internally in several parts of the code when parsing
things like numbers which are parts of strings (i.e. the 1 in
\inlinecode{direction1} for distribution functions). It probably shouldn't
be used to directly parse input deck parameters, since there is no reason to
restrict the user from specifying mathematical expressions.

\pagebreak
\begin{codedef}
FUNCTION as_long_integer_simple(str_in, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
INTEGER(KIND=8) :: as_long_integer_simple
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{as\_long\_integer\_simple} is equivalent to
\inlinecode{as\_integer\_simple}, but returns the larger
\inlinecode{INTEGER(KIND=8)} rather than a normal \inlinecode{INTEGER(KIND=4)}.
\\[0.5cm]
{\Large Notes\\}

\pagebreak
\begin{codedef}
FUNCTION as_boundary(str_in, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
INTEGER :: as_boundary
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{as\_direction} is used when assigning a laser to a boundary and
recognises the strings
\begin{itemize}
\item x\_min or left - c\_bd\_x\_min.
\item x\_max or right - c\_bd\_x\_max.
\item y\_min or down - c\_bd\_y\_min.
\item y\_max or up - c\_bd\_y\_max.
\item z\_min or back - c\_bd\_z\_min.
\item z\_max or front - c\_bd\_z\_max.
\end{itemize}
It returns the associated direction code (given after the dash in the
definition).
\\[0.5cm]
{\Large Notes\\}
If you're writing code which requires attaching something to a boundary,
whether a boundary condition, a diagnostic or some other routine, then this is
the routine that should be used. Note that in order to prevent confusion when
moving input decks between different dimension versions of {\EPOCH}, each code
only recognises the strings for boundaries that it actually has.

\pagebreak
\begin{codedef}
FUNCTION as_logical(str_in, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
Logical :: as_logical
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{as\_logical} simply tests for the strings ``T'' and ``F'' to
determine a boolean value. The default behaviour of \inlinecode{as\_logical} is
to treat any string that isn't ``T'' as a false value.
\\[0.5cm]
{\Large Notes\\}
You should use this rather than using a 0/1 boolean flag in the input deck for
consistency.

\pagebreak
\begin{codedef}
SUBROUTINE split_off_int(str_in, str_out, int_out, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
CHARACTER(LEN=*), INTENT(OUT) :: str_out
INTEGER, INTENT(OUT) :: int_out
INTEGER, INTENT(INOUT) :: err
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings\_advanced.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{split\_off\_int} is a routine which splits a string of the format
{\bf string{\it n}} into a string {\bf string} and an integer {\it n} which are
returned separately in the \inlinecode{str\_out} and \inlinecode{int\_out}
parameters respectively.  If it can't split the string successfully then it
sets the \inlinecode{c\_err\_bad\_value} bitfield of the err parameter.
\\[0.5cm]
{\Large Notes\\}
This is used in the core of the deck parser to deal with blocks like the
numbered species blocks in the initial conditions, and also in some of the
specific block parsers. Again, this routine should be used to split strings
like this rather than coding a new routine.

\pagebreak
\begin{codedef}
SUBROUTINE split_range(str_in, real1, real2, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
REAL(num), INTENT(OUT) :: real1, real2
INTEGER, INTENT(INOUT) :: err
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings\_advanced.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{split\_range} is a routine which splits a string of the format
{(\bf{\it n}, {\it m})} into two reals {\it n} and {\it m} which are returned
separately in the \inlinecode{real1} and \inlinecode{real2} parameters
respectively.  If it can't split the string successfully then it sets the
\inlinecode{c\_err\_bad\_value} bitfield of the err parameter.
\\[0.5cm]
{\Large Notes\\}
This is used when specifying ranges in the input decks at present. Any
ranges which should be specified in a single parameter should be specified in
this form and this routine used to split the string.\\ {\bf IMPORTANT NOTE. This
routine has proved somewhat unsatisfactory and may be replaced by another
routine in the future. It is expected that the function will have the same name
and behaviour, but might have subtle differences.}

\pagebreak
\begin{codedef}
FUNCTION as_integer(str_in, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
INTEGER :: as_integer
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings\_advanced.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{as\_integer} is the routine which returns integers from strings
using the maths parser. If a mathematical expression resolves to a non-integer
result then this routine rounds to the NEAREST integer. There are explicit
rounding routines in the maths parser to force other behaviour.
\\[0.5cm]
{\Large Notes\\}
This routine should be used when evaluating most strings into integers. Note
that this routine evaluates spatially dependent quantities at (0,0) on each
processor, so will give unpredictable results when spatially dependent
quantities are given to it (like density, bx etc.). To evaluate a spatially
varying quantity use \inlinecode{evaluate\_string\_in\_space}.

\pagebreak
\begin{codedef}
FUNCTION as_long_integer(str_in, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
INTEGER(KIND=8) :: as_long_integer
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings\_advanced.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{as\_long\_integer} is the routine which returns long integers from
strings using the maths parser. If a mathematical expression resolves to a
non-integer result then this routine rounds to the NEAREST integer. There are
explicit rounding routines in the maths parser to force other behaviour.
\\[0.5cm]
{\Large Notes\\}
This routine should be used when evaluating strings which are likely to be too
large to be stored in an INTEGER(KIND=4).

\pagebreak
\begin{codedef}
FUNCTION as_real(str_in, err)
\HRule
CHARACTER(LEN=*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
REAL(num) :: as_real
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings\_advanced.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{as\_real} is the routine which returns reals from strings using the
maths parser.
\\[0.5cm]
{\Large Notes\\}
This routine should be used when evaluating most strings into reals. Note that
this routine evaluates spatially dependent quantities at (0,0) on each
processor, so will give unpredictable results when spatially dependent
quantities are given to it (like density, bx etc.). To evaluate a spatially
varying quantity use \inlinecode{evaluate\_string\_in\_space}.

\pagebreak
\begin{codedef}
SUBROUTINE evaluate_string_in_space(str_in, data_out, &
    x1, x2, \{y1, y2, z1, z2,\} err)
\HRule
CHARACTER(*), INTENT(IN) :: str_in
INTEGER, INTENT(INOUT) :: err
INTEGER, INTENT(IN) :: x1, x2\{, y1, y2, z1, z2\}
REAL(num), DIMENSION(1:,1:,1:), INTENT(OUT) :: data_out
\end{codedef}
\vspace{1cm}
{\Large Location\\}
\inlinecode{src/deck/strings\_advanced.f90}\\
\\[0.5cm]
{\Large Description\\}
\inlinecode{evaluate\_string\_in\_space} is a routine which is used to evaluate
a tokenized maths expression over a region of the domain. The dimensionality of
\inlinecode{data\_out}, and the presence or absence of \inlinecode{y1, y2} and
\inlinecode{z1, z2} depend on the dimensionality of the code being used. The
\inlinecode{x1, x2, ...} parameters represent the indices in that
direction over which the expression should be evaluated. For example, in 2D to
evaluate an expression over the entire domain, the code would look like:
\begin{boxverbatim}
REAL(num), DIMENSION(:,:), ALLOCATABLE :: data
ALLOCATE(data(-2:nx+3,-2:ny+3)
CALL evaluate_string_in_space(string, data, -2, nx+3, -2, ny+3, err)
\end{boxverbatim}
{\Large Notes\\}
This routine is suitable to evaluate expressions over a subsection or all of
the domain, and is used in this way in the initial condition deck parser
routines. However, the routine does have one significant weakness, which is
that it tokenizes the string each time it is called. Tokenizing the string is a
time consuming process, so if the string is to be evaluated several times for
different reasons (for example, the time profile for the laser) then a
different procedure should be followed using the lower level parser routines,
which are described in \sect{maths}.
\pagebreak

\subsection{Adding new elements to the maths parser}
Sometimes the complexity in changing the input.deck file is due to the fact that
a function which must be used is fairly complex in form and is not supplied
with the core code. It must therefore be represented in the input deck maths
parser. This can be a significant cause of complexity for some problems, and
in this case, there are three options: Put up with it and implement in the
deck, use the internal initial conditions rather than the deck or extend the
maths parser to include your function. Extending the maths parser can either
be permanent (described later in \sect{maths}) or temporary (described
here). Temporarily adding elements to the parser is much the easier. It is
possible to add new constants and functions to the maths parser. It is hoped
that in a future release of {\EPOCH} this will be extended to allow custom
operators as well.\\

As an example, lets look at adding a new function (lorentz) for a
Lorentzian distribution, and adding a new constant, phi.

\subsubsection{Registering your new constant/function}
All of the routines used in extending the maths parser are in the file
user\_interaction/custom\_parser.f90.  Before a new constant or function
can be defined it must be registered. In the registration phase the text
representation of the function or constant is given to the parser subroutines
and the user is returned an integer handle for the registered object. The
numerical handle must be stored so that that all of the functions in this
module can access it, so they should be placed after the \inlinecode{IMPLICIT
NONE} statement at the top of the file and defined as:
\begin{boxverbatim}
INTEGER :: c_func_lorentz
INTEGER :: c_const_phi
\end{boxverbatim}

Note that the names given to the constants is obviously at the developers
discretion, but these names comply with the {\EPOCH} style guide.
Actually registering the objects is done in the \inlinecode{register\_objects}
subroutine which should include lines to register functions and constants.
An example is given below.
\begin{boxverbatim}
SUBROUTINE register_objects

  c_func_lorentz = register_function("lorentz")
  c_const_phi = register_constant("phi")

END SUBROUTINE register_objects
\end{boxverbatim}

Note that the input deck parser is case sensitive, so the strings which are
given to \inlinecode{register\_function} and \inlinecode{register\_constant}
should be in the case that they will appear in the input deck. To follow the
{\EPOCH} style guide this should be all lowercase. At this point, the maths
parser would start to recognise the new function/constant, but would still
give error messages since they haven't yet been implemented.

\subsubsection{Setting up new constants}

Once a new constant has been registered it must be described using the
\inlinecode{custom\_constant} function. In 2D this function looks like:
\begin{boxverbatim}
FUNCTION custom_constant(opcode, ix, iy, errcode)

  INTEGER, INTENT(IN) :: opcode, ix, iy
  INTEGER, INTENT(INOUT) :: errcode
  REAL(num) :: custom_constant

  ! Leave these lines in place. They cause the code to throw an error if
  ! The opcode is unknown
  custom_constant = 0.0_num
  errcode = IOR(errcode, c_err_unknown_element)

END FUNCTION custom_constant
\end{boxverbatim}

The parameters are

\begin{itemize}
\item opcode - The operator code of the constant requested. This will be the
  integer handle returned from \inlinecode{register\_constant}.
\item ix, iy, iz - Some constants are actually evaluated at specific points in
  space and ix, iy, iz are the gridpoint number of the location currently being
  evaluated. If you are specifying a simple constant then just ignore
  these. If your constant does depend upon space then directly subscript your
  array with ix, iy, iz as needed to read the correct location.
\item errcode - The error code which should be passed back to the
  parser. If for some reason you cannot evaluate your constant then you should
  \inlinecode{IOR} errcode with the appropriate error code (all the error
  codes are listed in appendix A). Note that errcode should never be SET to
  any specific error code when extending the parser, since this might
  overwrite errors put in place earlier in the parsing sequence. This is
  different to extending the input deck where the error code is set.
\end{itemize}

The function should just return the evaluated value of the constant requested
by \inlinecode{opcode}. This might look like:
\begin{boxverbatim}
FUNCTION custom_constant(opcode, ix, iy, errcode)

  INTEGER, INTENT(IN) :: opcode, ix, iy
  INTEGER, INTENT(INOUT) :: errcode
  REAL(num) :: custom_constant

  IF (opcode .EQ. c_const_phi) THEN
    custom_constant = pi
    RETURN
  ENDIF

  ! Leave these lines in place. They cause the code to throw an error if
  ! The opcode is unknown
  custom_constant = 0.0_num
  errcode = IOR(errcode, c_err_unknown_element)

END FUNCTION custom_constant
\end{boxverbatim}

Note that when \inlinecode{opcode} is successfully recognised, the code sets
the return value and returns straight away. This is how all constants should
work, since the last line forces the function to return an error code. This
last line is in place to trap people registering constants but never defining
them. Without this line, it is possible to define a constant which is
never specified and have the code complete OK with a random value for that
constant.

The constant ``phi'' should now work fine when used anywhere in the input deck
and will return a value of $\pi$.

\subsubsection{Setting up new functions}

Setting up the new function \inlinecode{lorentz} is very similar to setting up
the new constant. The relevant function is \inlinecode{custom\_function} and
when empty looks like:
\begin{boxverbatim}
FUNCTION custom_function(opcode, ix, iy, errcode)

  INTEGER, INTENT(IN) :: opcode, ix, iy
  INTEGER, INTENT(INOUT) :: errcode
  REAL(num) :: custom_function
  REAL(num) :: values(5)

  ! Leave these lines in place. They cause the code to throw an error if
  ! The opcode is unknown
  custom_function = 0.0_num
  errcode = IOR(errcode, c_err_unknown_element)

END FUNCTION custom_function
\end{boxverbatim}

The parameters are

\begin{itemize}
\item opcode - The operator code of the constant requested. This will be the
  integer handle returned from \inlinecode {register\_function}.
\item ix, iy, iz - Some functions are evaluated differently at specific points
  in space and ix, iy, iz are the gridpoint number of the location currently
  being evaluated. If you are specifying a simple function then just
  ignore these. If your function does depend upon space then directly
  subscript your array with ix, iy, iz as needed to read the correct location.
\item errcode - The error code which should be passed back to the
  parser. If for some reason you cannot evaluate your function then you should
  \inlinecode{IOR} errcode with the appropriate error code. Note that errcode
  should never be SET to any specific error code, since this might overwrite
  errors put in place earlier in the parsing sequence.
\end{itemize}

The function should return the value of your evaluated constant. The
parameters which are passed to the function can be retrieved by the function
\inlinecode {get\_values(n, values)}, where \inlinecode {n} is the number of
parameters to be returned and \inlinecode{values} is a \inlinecode{REAL(num)}
array of length \inlinecode{n} which will hold the returned values .  In this
implementation of the Lorentzian function there are three parameters: The
dependent variable, the location parameter and the scale parameter. The code
to implement the function therefore looks like:
\begin{boxverbatim}
FUNCTION custom_function(opcode, ix, iy, errcode)

  INTEGER, INTENT(IN) :: opcode, ix, iy
  INTEGER, INTENT(INOUT) :: errcode
  REAL(num) :: custom_function
  REAL(num) :: values(5)

  IF (opcode .EQ. c_func_lorentz) THEN
    CALL get_values(3, values(1:3))
    ! values(1) - Dependent variable
    ! values(2) - location parameter
    ! values(3) - scale parameter
    custom_function = values(3)**2/((values(1)-values(2))**2+values(3)**2)
    RETURN
  ENDIF

  ! Leave these lines in place. They cause the code to throw an error if
  ! The opcode is unknown
  custom_function = 0.0_num
  errcode = IOR(errcode, c_err_unknown_element)

END FUNCTION custom_function
\end{boxverbatim}

This function is then available at any point in the input deck and if I return
to the previous example ic.deck file, it would be used as follows:
\begin{boxverbatim}
begin:constant
   particle_density = 100.0 # Particle number density
   profile_x = lorentz(x,0.0,1.0)
   profile_y = lorentz(y,0.0,1.0)
end:constant

begin:species
   name = s1

   # multiply density by real particle density
   density = particle_density * profile_x * profile_y

   # Set the temperature to be zero
   temp_x = 0.0
   temp_y = temp_x(s1)

   # Set the minimum density for this species
   density_min = 0.3*particle_density
end:species

begin:species
   name = s2

   # Just copy the density for species s1
   density = density(s1)

   # Just copy the temperature from species s1
   temp_x = temp_x(s1)
   temp_y = temp_y(s1)

   # Set the minimum density for this species
   density_min = 0.3*particle_density
end:species
\end{boxverbatim}

It is therefore clear that the new lorentz function is essentially the same as
the built in gauss function. Note that due to the way that the parser works,
the end user is not required to deal with parameters which are themselves
maths expressions. They have been fully evaluated by the time they are
returned by \inlinecode{get\_values}. Note that the parser is not guaranteed to
be bulletproof. If a user calls \inlinecode{get\_values} requesting more
parameters than have been passed to the function then it will scramble the
stack which is used by the parser and cause the code to fail. Note that
calling \inlinecode{get\_values(2, values)} is not the same as calling
\inlinecode{get\_values(1, values)} twice, in fact calling
\inlinecode{get\_values(1, values)} multiple times will return the parameters in
{\it reverse} order. This is normal and is a feature of how the maths parser
operates. It is possible to use this property to write functions which have a
variable number of parameters, but this is not recommended.

\subsection{Adding new elements to the input deck}

For some types of changes to the code it is more convenient to have the end
user pass new parameters into the code. This can be for several reasons, and
the section on permanent additions to the input deck is given later in this
manual. At this stage, we will describe how to temporarily add new
elements to the input deck parser routines, allowing parameterising of
internal and manual initial conditions.

Custom input deck elements are setup in the file
\inlinecode{src/user\_interaction/custom\_deck.f90}. The function
\inlinecode{custom\_blocks\_handle\_element} is called when a new block is
started which the core parser is not familiar with, and once for each element
of a block. The function \inlinecode{custom\_blocks\_check} is called once the
entire deck has been parsed and is used to check that all the elements which
are required for the code to run have been specified.

\subsubsection{custom\_blocks\_handle\_element}
There are three parameters passed to
\inlinecode{custom\_blocks\_handle\_element}, which are:
\begin{itemize}
\item block\_name - The name of the block specified in the
  \inlinecode{begin:blockname} part of the input deck.
\item element - The name of the element in an input deck
  \inlinecode{element = value} pair.
\item value - The string representation of the value in an input deck
  \inlinecode{element = value} pair.
\end{itemize}

\inlinecode{custom\_blocks\_handle\_element} is first called when a new block
is begun
using \inlinecode{begin:blockname} and ``blockname'' is not recognised by the
core input deck parser. The first thing that it does is test whether or not it
is a valid custom block. The code does this by passing in the blockname with
\inlinecode{element} and \inlinecode{value} set to the special constant called
``blank''. When extending the input deck, an end user should check if either
\inlinecode{element} or \inlinecode{value} are set to the special constant
``blank'', and if they are then test to see whether the blockname is known or
not. If the blockname is known then the code should return the error code
\inlinecode{c\_err\_none} (No error). If the blockname is not known then the
code should return \inlinecode{c\_err\_unknown\_block} and the deck parser will
just skip the block. In operation, this looks like:
\begin{boxverbatim}
FUNCTION custom_blocks_handle_element(block_name, element, value) &
    RESULT(errcode)

  CHARACTER(LEN=string_length), INTENT(IN) :: block_name, element, value
  INTEGER :: errcode

  IF (str_cmp(block_name, "custom")) THEN
    IF (element .EQ. blank .OR. value .EQ. blank) THEN
       ! If element or value are blank then just testing block so
       ! return c_err_none
       errcode = c_err_none
       RETURN
     ENDIF
   ENDIF

  ! The following line must always be present
  errcode = c_err_unknown_block

END FUNCTION custom_blocks_handle_element
\end{boxverbatim}

In order to simplify Fortran's rather annoying string handling behaviour,
several helper functions have been defined and the most used one is
\inlinecode{str\_cmp(string1, string2)}. This is a simple routine which returns
true if string1 == string2 and false otherwise. It is case sensitive but can
deal with differing string lengths etc. The next stage is to deal with the
actual \inlinecode{element = value} pairs in the deck. Each time that a new pair
is read from the deck, \inlinecode{custom\_blocks\_handle\_element} is called
with \inlinecode{element} and \inlinecode{value} having the values read from the
deck. To test for known elements they should just be checked against a known
list of names using \inlinecode{str\_cmp} and return the error code
\inlinecode{c\_err\_unknown\_element} if the element isn't a known element. This
looks like:
\begin{boxverbatim}
FUNCTION custom_blocks_handle_element(block_name, element, value) &
    RESULT(errcode)

  CHARACTER(LEN=string_length), INTENT(IN) :: block_name, element, value
  INTEGER :: errcode

  IF (str_cmp(block_name, "custom")) THEN
    IF (element .EQ. blank .OR. value .EQ. blank) THEN
      ! If element or value are blank then just testing block so
      ! return c_err_none
      errcode = c_err_none
      RETURN
    ENDIF
    errcode = c_err_unknown_element
    ! Now test for the real elements
    IF (str_cmp(element, "int_element")) THEN
      errcode = c_err_none
    ENDIF
    IF (str_cmp(element, "real_element")) THEN
      errcode = c_err_none
    ENDIF
    IF (str_cmp(element, "logical_element")) THEN
      errcode = c_err_none
    ENDIF
    RETURN
  ENDIF

  ! The following line must always be present
  errcode = c_err_unknown_block

END FUNCTION custom_blocks_handle_element
\end{boxverbatim}

This version of the code will allow you to add a new block called ``custom''
with elements
``int\_element'', ``real\_element'' and ``logical\_element'' and the
code will parse them successfully, while any other block or any other element
in the block ``custom'' will throw errors. However, at this stage the code
doesn't actually read any of the values from the deck. To make it useful, any
variable which is read from the input deck must be stored in a global
variable. Defining global variables are explained in more detail in the
relevant section of the manual, but in short, any variable defined in the
module \inlinecode{shared\_data} in the file \inlinecode{src/shared\_data.F90}
will be a global variable. After the variables have been setup, there are once
again helper functions to make converting the text from the deck into a normal
Fortran90 variable. These helper functions are:

\begin{itemize}
\item as\_integer - Attempts to convert a string to an integer. Invokes the
  maths parser.
\item as\_real - Attempts to convert a string to a REAL(num). Invokes the maths
  parser.
\item as\_logical - Attempts to convert a string to a logical. Does not invoke
  the maths parser (must be either ``T'' or ``F'').
\end{itemize}

They are used pretty much as expected, except that the return value is passed
to the functions so that they can report errors while trying to parse the
string. An example would then be:
\begin{boxverbatim}
FUNCTION custom_blocks_handle_element(block_name, element, value) &
    RESULT(errcode)

  CHARACTER(LEN=string_length), INTENT(IN) :: block_name, element, value
  INTEGER :: errcode, int_element
  REAL(num) :: real_element
  LOGICAL :: logical_element

  IF (str_cmp(block_name, "custom")) THEN
    IF (element .EQ. blank .OR. value .EQ. blank) THEN
      ! If element or value are blank then just testing block so
      ! return c_err_none
      errcode = c_err_none
      RETURN
    ENDIF
    errcode = c_err_unknown_element
    ! Now test for the real elements
    IF (str_cmp(element, "int_element")) THEN
      errcode = c_err_none
      int_element = as_integer(value, errcode)
    ENDIF
    IF (str_cmp(element, "real_element")) THEN
      errcode = c_err_none
      real_element = as_real(value, errcode)
    ENDIF
    IF (str_cmp(element, "logical_element")) THEN
      errcode = c_err_none
      logical_element = as_logical(value, errcode)
    ENDIF
    RETURN
  ENDIF

  ! The following line must always be present
  errcode = c_err_unknown_block

END FUNCTION custom_blocks_handle_element
\end{boxverbatim}

It is possible to perform more advanced types of evaluation of maths
expressions such as reading arrays etc. but this is beyond the scope of this
manual at present.

\subsubsection{custom\_blocks\_check}
This function is called when all the blocks in the input deck have been
evaluated and is used to check that all required parameters have been set. If
all required elements have been set then you should just return
\inlinecode{c\_err\_none}, otherwise return
\inlinecode{c\_err\_missing\_elements}. How you test that required elements have
been set is up to the developer, and for testing and personal use (which is
all that the custom deck parts of the code should be used for) it is
acceptable to just not check and always return \inlinecode{c\_err\_none}. If
permanently expanding the deck, error trapping should always be written.

\subsection{Permanently adding blocks to the input deck}

While using the \inlinecode{custom\_deck} subroutine is a good way of passing
parameters into the code or for temporary additions, it is not suitable for
permanent additions to the code. Adding new blocks to the code permanently is
very similar to doing it temporarily, but requires changes to some of the
subroutines in \inlinecode{deck.F90}.\\

There are six subroutines which may need to be changed to add new blocks to
the deck. These are
\begin{itemize}
\item \inlinecode{deck\_initialise} - Called before parsing of the input
  deck is begun.
\item \inlinecode{deck\_finalise} - Called after parsing of the input
  deck is complete.
\item \inlinecode{start\_block(block\_name)} - Called when the deck directive
  \inlinecode{begin:{\it block\_name}} appears in a deck file.
\item \inlinecode{end\_block(block\_name)} - Called when the deck directive
  \inlinecode{end:{\it block\_name}} appears in a deck file.
\item \inlinecode{handle\_block(block\_name, block\_element, block\_value)} -
  Called once for each element in a block.
\item \inlinecode{check\_compulsory\_blocks(errcode\_deck)} - Called once when
  the deck file has been read to check that all necessary blocks have been
  populated.
\end{itemize}

  %FUNCTION species_block_handle_element(element, value) RESULT(errcode)
  %FUNCTION species_block_check() RESULT(errcode)


There is one final variable which is important for modifying the input deck,
\inlinecode{deck\_state}. The input deck parser routine used to read the main
input deck uses the variable \inlinecode{deck\_state} to
determine which stage of parsing the deck is required. The possible values of
\inlinecode{deck\_state} are

\begin{itemize}
\item c\_ds\_first - The first pass through the deck, before memory has been
  allocated.
\item c\_ds\_last - After the initial deck pass, all arrays and lists are
  allocated. The deck is then parsed a final time so that allocated memory
  can be populated with initial conditions.
\end{itemize}
These constants are defined in \inlinecode{shared\_data.F90}.

The layout of \inlinecode{deck\_initialise} and \inlinecode{deck\_finalise} is
extremely simple. They just call \inlinecode{*\_deck\_initialise} or
\inlinecode{*\_deck\_finalise} for each of the possible block types.
\inlinecode{start\_block} and \inlinecode{end\_block} are also fairly
straightforward. They examine the block name and call the
\inlinecode{*\_block\_start} or \inlinecode{*\_block\_end} routine
corresponding to the current block.

The \inlinecode{handle\_block} routine acts in a similar manner except
that it also does some error handling.
At the simplest level the routine simply calls another function which
takes the block\_element and block\_value as
parameters and returns an error code from Appendix~\ref{sec:error}
determining the success or failure of reading the element.

The final routine is \inlinecode{check\_compulsory\_blocks} which is used to
check that all the needed elements of the input deck have been set. A single
parameter \inlinecode{errcode\_deck} is passed in. The routine
checks \inlinecode{deck\_state} to make sure that it is the last pass
through the deck. It then goes through and calls functions to check that
all the necessary parts of a block have been set. The subroutines are contained
in the same file as the routine which is called in \inlinecode{handle\_block} to
handle elements of the block. The error handler functions should return an
error code from Appendix~\ref{sec:error}, usually
\inlinecode{c\_err\_missing\_elements}. The
return code from the error handler function should then be \inlinecode{IOR}ed
with \inlinecode{errcode\_deck} to allow error codes to be returned from
several different checks with errors occurring.

\subsubsection{The element handler routines for deck elements}
The exact form of the handler routines is up to the end user. The only {\it
requirements} are that the routine should return an error code detailing
whether or not there are any problems with reading the block and that the error
code should be \inlinecode{c\_err\_none} if either the element name or element
value are the special constant \inlinecode{blank}. The typical implementation
of an element handler routine is shown in the file
\inlinecode{src/deck/deck\_control\_block.f90}, and this general layout should
be copied for compatibility if possible.\\

Sometimes, it is useful to have each new block correspond to a new instance of
an object in the code. An example of this in {\EPOCH} is in
\inlinecode{src/deck/deck\_laser\_block.f90} where each new laser block in
the input deck corresponds to a new laser being attached to a boundary. This is
accomplished by implementing the lasers as a linked list on each boundary,
with a new laser object being created when a laser block is started, the laser
information being set during the main reader routine, and then the laser being
attached to the linked list by a call to \inlinecode{attach\_laser} in
\inlinecode{src/laser.f90} when the block is ended. When a new laser block is
started the process simply repeats allowing the end user to have as many lasers
as desired.

\subsubsection{Adding elements to existing blocks}
The existing blocks in the code are read in the files listed in \sect{src_deck}

The existing structure of the blocks is simple enough in most cases that it
should be fairly easy to add new elements if needed. The most likely change
needed is to change the list of variables to dump in the \inlinecode{output}
block. How to do this is detailed in \sect{io}.

\subsection{Permanently adding functions, constants to {\EPOCH}}
\label{sec:maths}
In much the same way that it is possible to permanently add new blocks to the
input deck, it is also possible to permanently add functions and constants to
{\EPOCH}'s maths parser. Although adding new operators is possible, it is
sufficiently likely to cause problems with the operation of the maths parser
that it is formally not recommended by the author of the program, and hence is
not documented here.

\subsubsection{Adding the new tokenizer handle}
When adding a new function or constant to the maths parser using the temporary
routines, there are two calls (\inlinecode{register\_function} and
\inlinecode{register\_constant}) which give a numerical handle. This is the
token used to represent that function or constant after the text has been parsed
(remember that {\EPOCH}'s maths parser tokenizes before evaluation!). When
permanently adding objects to the maths parser, the tokenizer handles have to
be set up manually. This takes place in \inlinecode{src/shared\_data.F90} in
the module \inlinecode{shared\_parser\_data}. There are several lines which
look like:
\begin{boxverbatim}
  INTEGER, PARAMETER :: c_const_ix = 40
  INTEGER, PARAMETER :: c_const_iy = 41
  .
  .
  .
  INTEGER, PARAMETER :: c_func_interpolate = 22
  INTEGER, PARAMETER :: c_func_tanh = 23
\end{boxverbatim}
Constants beginning with \inlinecode{c\_const\_} are tokenizer handles for
constants, and those beginning with \inlinecode{c\_func\_} are tokenizer handles
for functions. Each number must be unique and has to be less than
the lower bound of values reserved for temporary or deck
specified values. This means that any tokenizer handle for a function has to be
less than the value of the variable \inlinecode{c\_func\_custom\_lowbound} and
any handle for a constant must be less than
\inlinecode{c\_const\_deck\_lowbound}. It is acceptable to simply increase the
value of \inlinecode{c\_func\_custom\_lowbound} and
\inlinecode{c\_const\_deck\_lowbound} to
allow the use of more values for internal constants and functions, although
care should be taken since this could lead to performance issues.
If \inlinecode{c\_const\_deck\_lowbound} is increased then the constant
\inlinecode{c\_constant\_custom\_lowbound} should be increased by the same
amount (the values between \inlinecode{c\_const\_deck\_lowbound} and
\inlinecode{c\_constant\_custom\_lowbound-1} are used for constants specified
inside the input deck while values greater than or equal to
\inlinecode{c\_constant\_custom\_lowbound} are used for constants specified
by \inlinecode{register\_constant}.

Once the tokenizer handle is specified in \inlinecode{shared\_parser\_data}, it
is now possible to extend the main areas of the maths parser. Note that from
here on, you MUST always use the constant named handle, NEVER the numerical
value that you specified for the value of the handle. If this is not done
then combining functions and constants from several sources becomes much harder.

\subsubsection{Adding the new function or constant to the tokenizer}
The next stage is to add the string representation of your constant or function
to the tokenizer routines in
\inlinecode{src/parser/tokenizer\_blocks.f90}. This is very simple to do, just
find either the function \inlinecode{as\_constant} or \inlinecode{as\_function}
and look at the existing code. These functions are just long lists of
\inlinecode{str\_cmp} commands followed by code to deal with custom
functions/constants. To add the new code, create an additional line such as:
\begin{boxverbatim}
  IF (str_cmp(name, "my_const")) as_constant = c_const_my_const
  .
  .
  .
  IF (str_cmp(name, "my_func"))  as_function = c_func_my_func
\end{boxverbatim}
Note that neither routine returns immediately after recognising the name of the
function/constant. This allows users to override built in constants or
functions with custom versions using \inlinecode{register\_constant}
and \inlinecode{register\_function}. This is not significant, since tokenizing
should never be used in a speed critical part of the code.

\subsubsection{Implementing the function or constant in the evaluator}
The evaluator is the part of the code that actually takes the streams of tokens
produced by the tokenizer and evaluates them into a number. The relevant parts
of the evaluator for adding new constants or functions are in
\inlinecode{src/parser/evaluator\_blocks.f90} and the functions which may need
changing are \inlinecode{do\_constant} and \inlinecode{do\_functions}. These are
both passed up to five parameters:
\begin{itemize}
\item #INTEGER :: opcode# - The operation code, this is the tokenizer handle
  which was defined in \inlinecode{shared\_parser\_data}.
\item #INTEGER :: ix, iy, iz# - The position of the current evaluation in the
  domain. If your function or constant behaves differently at different points
  in space then you should use these parameters to reference the correct point
  of an array.
\item #INTEGER :: errcode# - This should be set to an error code from
  Appendix~\ref{sec:error}, usually \inlinecode{c\_err\_bad\_value} if for some
  reason it is not possible to evaluate your constant or function.
\end{itemize}

The rest of the routine to set a constant is as simple as testing for the
tokenizer handle already set up in \inlinecode{shared\_parser\_data} and then
calling the subroutine \inlinecode{push\_on\_eval}. This pushes the final
constant onto the evaluation stack which is used by the RPN parser. The basic
sequence for functions is similar except for the addition of a code to read
the values that the function takes. This is again the subroutine
\inlinecode{get\_values} which is also used in custom\_function. The calling
sequence in \inlinecode{do\_function} looks like:
\begin{boxverbatim}
  IF (opcode .EQ. c_func_gauss) THEN
    CALL get_values(3, values)
    CALL push_on_eval(EXP(-((values(1)-values(2))/values(3))**2))
    RETURN
  ENDIF
\end{boxverbatim}
Simply call the \inlinecode{get\_values} subroutine, passing the number of
required parameters and an array of type \inlinecode{REAL(num)} which is at
least as long as the number of required parameters. The array is populated
by the values passed into the function. Constants and maths expressions are
already evaluated by the time that this section of code is reached, so there is
no need to deal with further parsing. Next, simply call
\inlinecode{push\_on\_eval} to push the result of your function onto the
evaluation stack.

\section{Developing an extension to {\EPOCH}}
\label{sec:extend}

Exactly how to extend {\EPOCH} depends heavily upon what you intend to add. The
simplest things to add are new diagnostics, and this is detailed in
\sect{io}. Other places where changes are likely to be made
are the following:
\begin{itemize}
\item The field solver - Changing the field solver to add new laser-like
  boundaries, add spatial smoothing to remove noise, add high order field
  solvers etc.
\item The particle pusher - Change the basic physical model of {\EPOCH} by
  modifying the particle pusher.
\item The boundary routines - Add new boundary conditions or modify existing
  boundary conditions.
\item The laser boundary routines - Add new features to the laser boundaries in
  this routine.
\item The main driver (\inlinecode{epoch\{n\}d.F90}) - This is the routine
  where the main calling sequence of {\EPOCH} is setup, and totally new
  extensions to {\EPOCH} should be placed in here.
\end{itemize}

Changing the field solver, the particle pusher or boundary routines
is fairly easy to accomplish by reading the section of this manual that
details the relevant code. The general sequence for writing an extension
would be:
\begin{itemize}
\item Add any new global variables needed to \inlinecode{shared\_data.F90}.
\item Add the meat of your change to the code.
\item Test the changes to your code. Make absolutely sure that you can turn your
  change to the code off.
\item Add controls for your extension to the input deck reader.
\end{itemize}

\subsection{The main driver routine}
When adding completely new routines to the code, they should be added to the
file \inlinecode{src/epoch\{n\}d.F90}. This routine simply calls other routines
to perform the actual execution of the code. The first section of the code
controls basic setup, MPI initialisation and
initial conditions. If you wish to add new startup conditions then
you should find the location in this routine where the
initial conditions are setup. The code is fairly complicated, but
there are a few key points at which the code significantly changes state.
\begin{itemize}
\item After the call to \inlinecode{read\_deck} the code has read the basic
  information from the input deck files and any tests or changes which have to
  be made to input deck values should be made immediately after this line. Note
  that although the variables from the deck have been set,
  none of these values have been used so allocatable
  variables have yet to be allocated. The grid does not exist at this point.
\item After the call to \inlinecode{open\_files} the code has finished
  allocating all field variables, although particles may not yet have been set
  up. The grid now exists.
\item There are now a series of \inlinecode{IF} statements which test for
  things like \inlinecode{IF (IOR(ictype, c\_ic\_autoload) .NE. 0)}.
  These are the lines which test for all possible states of the
  initial conditions. The last test is for the manual load routine
  (\inlinecode{c\_ic\_manual}). After this test all the particles have been
  loaded and are now on their correct processor. The load balancer has now been
  called at least once so the domains may no longer be identical.
\item The main loop is a simple do loop beginning with just the single command
  \inlinecode{DO}. Inside this loop there are several calls to routines which
  actually advance the system. Most routines which can change currents should
  take place after the particle pusher but before the final update for the $E$
  and $B$ fields. These routines are
  \begin{itemize}
  \item \inlinecode{set\_dt} - This routine sets the timestep.
  \item \inlinecode{update\_eb\_fields\_half} - Time centre the $E$ and $B$
    fields.
  \item \inlinecode{push\_particles} - The particle pusher.
  \item \inlinecode{reorder\_particles\_to\_grid} - Groups particles into
    linked lists at each grid point. Used for the particle splitting routine,
    and would be the right place to add a collision operator. Any routine
    which needs to have nearby particles grouped together should take place
    after the call to this routine.
  \item \inlinecode{split\_particles} - The particle
    splitting operator.
  \item \inlinecode{reattach\_particles\_to\_mainlist} - Undoes the
    particle grouping and rebuilds the main list of particles used by the
    particle pusher. Any routine which needs to have nearby particles grouped
    together should take place before the call to this routine.
  \item \inlinecode{update\_eb\_fields\_final} - Updates the $E$ and $B$
    fields to the full timestep.
  \end{itemize}
\item After the call to \inlinecode{update\_eb\_fields\_final} the code is
  ready for another timestep. Any routines which do not change the time
  integrated properties of the code (like the moving window) should come after
  this call.
\end{itemize}

\subsection{The particle reordering routine}
If the code is compiled with the right flags then during the main loop the
particles are grouped into a linked list for each cell in the domain, with all
the particles which are in that cell linked into the list. The main list
\inlinecode{species(ispecies)\%attached\_list} is empty and cannot be used
during this period. The particles should now be accessed using the variable
\inlinecode{species(ispecies)\%secondary\_list(ix,iy,iz)} which is the array of
linked lists. This array is allocated on the call to
\inlinecode{reorder\_particles\_to\_grid} and deallocated on the call to
\inlinecode{reattach\_particles\_to\_mainlist}, and should not be used outside
the section of code between these two calls. The particles themselves remain
unchanged. No attempt is made to check that particles do not cross processor
boundaries in this section, so if a particle's position is changed, it is up to
the user to ensure that the particle is transferred to another processor if
required. However, if a particle is transferred to another processor, it is
acceptable to relink it to \inlinecode{species(ispecies)\%attached\_list} since
the other lists are simply appended to that list when the particles are
reattached to the main list.\\

You should not write an extension to the code which requires the use of the
particle reordering routines if it can be avoided since these routines
have a significant impact on performance.

\pagebreak
\def\epoch{}\input{sdf_format}
\pagebreak

\appendix
  \appendixpage
  \addappheadtotoc
\section{Error Codes}
\label{sec:error}
The input deck and maths parser in {\EPOCH} use various named error codes to
report on errors which occur during the evaluation of the input deck. These
codes are
\begin{itemize}
\item c\_err\_none - No error. Set an error code to c\_err\_none to state that
  no error has occurred.
\item c\_err\_unknown\_block - In the input deck a block has been found which is
  not known. This should be returned in
  \inlinecode{custom\_blocks\_handle\_element} if it is passed any block that it has not been written to handle.
\item c\_err\_unknown\_element - In the input deck an element of a valid block
  has been found which is not known. This should be returned in
  \inlinecode{custom\_blocks\_handle\_element} if an element is requested which is
  unknown.
\item c\_err\_preset\_element - An element of the input deck has already been
  set and is being set again. Usually this is an indication of a malformed input
  deck file, so \inlinecode{custom\_blocks\_handle\_element} should try to
  identify such situations and return this error message if
  subsequent attempts to set the variable are being ignored.
\item c\_err\_preset\_element\_use\_later - An element of the input deck has
  already been set and is being set again. Usually this is an indication of a
\setlength{\emergencystretch}{3em}
  malformed input deck file, so \inlinecode{custom\_blocks\_handle\_element} should try to
  identify such situations and return this error message if
  the subsequent attempts to set the variable override previous ones.
\item {\bf c\_err\_bad\_value} - A value which is being evaluated for the right
  hand side of an element assignment is in some way invalid. Internally to the
  code this usually means that a string which must be interpreted as a maths
  expression or numerical constant is in some way malformed. It is also
  acceptable to return this error code when a value has been passed which is
  invalid for some other reason (the value is outside an acceptable range, etc.)
\item {\bf c\_err\_missing\_elements} - This is an error code returned when the
  code is testing to make sure that all necessary elements of an input deck
  file have been specified. It should be returned when some required parameter
  is missing in the subroutine \inlinecode{custom\_blocks\_check}.
\item c\_err\_terminate - This error code means that the code is in a state
  where execution is impossible and the code must terminate once the input deck
  has been read. Some other error codes automatically set c\_err\_terminate,
  but it can always be IOR'ed with any error code to force the code to exit.
  Note that just returning c\_err\_terminate will cause the code to
  silently quit.
\item {\bf c\_err\_required\_element\_not\_set} - This means that the code
  cannot parse an input deck element since another element which must be known
  beforehand has not been set. This is intended for things like setting the
  species information where the number of species must be known in
  advance. This error code uses the extended error string to give user friendly
  feedback. If you return this error code then you should set
  extended\_error\_string to be equal to the name of the required element which
  has not been set. If multiple previous elements are required then the code
  should be set up so that it checks for the presence of the required elements
  in order and reports on missing elements so that the end user can fix them
  one by one.
\item c\_err\_pp\_options\_wrong - If you've written a section of the code that
  is controlled by preprocessor options then you should return this
  error message if someone attempts to set input deck elements which refer to
  that part when the correct preprocessor options are not used. This means that
  the user is aware of the fact that the requested feature will not be
  active. This error code also uses the extended error string to give user
  friendly feedback. If you return this error code, you should set the string
  extended\_error\_string to the define command that would turn on the
  requested feature of the code (``-DPER\_PARTICLE\_WEIGHT'', for example).
\item {\bf c\_err\_other} - This error code is a catch all error which causes
  the code to quit with a sarcastic error message. It's mainly intended for
  debugging and is used before the final error code is implemented.
\end{itemize}

Those error codes with names in {\bf BOLD} are those which will cause the code
to terminate.

\end{document}