{'fraction_non_alphanumeric': 0.038971576358426935, 'fraction_numerical': 0.010446512335505464, 'mean_word_length': 4.850763981605103, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Object-oriented implementations of the MPDATA advection equation solver in C++, Python and Fortran

Three object-oriented implementations of a prototype solver of the advection equation are introduced. The presented programs are based on Blitz++ (C++), NumPy (Python), and Fortran's built-in array containers. The solvers include an implementation of the Multidimensional Positive-Definite Advective Transport Algorithm (MPDATA). The introduced codes exemplify how the application of object-oriented programming (OOP) techniques allows to reproduce the mathematical notation used in the literature within the program code. A discussion on the tradeoffs of the programming language choice is presented. The main angles of comparison are code brevity and syntax clarity (and hence maintainability and auditability) as well as performance. In the case of Python, a significant performance gain is observed when switching from the standard interpreter (CPython) to the PyPy implementation of Python. Entire source code of all three implementations is embedded in the text and is licensed under the terms of the GNU GPL license.


Introduction
Object oriented programming (OOP) "has become recognised as the almost unique successful paradigm for creating complex software" [1,Sec. 1.3]. It is intriguing that, while the quoted statement comes from the very book subtitled The Art of Scientific Computing, hardly any (if not none) of the currently operational weather and climate prediction systems -flagship examples of complex scientific software -make extensive use of OOP techniques. Fortran has been the language of choice in oceanic [2], weather-prediction [3] and Earth system [4] modelling, and none of its 20-century editions were object-oriented languages [see e.g. 5, for discussion].
Application of OOP techniques in development of numerical modelling software may help to: (i) maintain modularity and separation of program logic layers (e.g. separation of numerical algorithms, parallelisation mechanisms, data input/output, error handling and the description of physical processes); and (ii) shorten and simplify the source code and improve its readability by reproducing within the program logic the mathematical notation used in the literature.
The first application is attainable, yet arguably cumbersome, with procedural programming. The latter, virtually impossible to obtain with procedural programming, is the focus of this paper. It also enables the compiler or library authors to relieve the user (i.e. scientific programmer) from hand-coding optimisations, a practice long recognised as having a strong negative impact when debugging and maintenance are considered [6].
MPDATA [7] stands for Multidimensional Positive Definite Advective Transport Algorithm and is an example of a numerical procedure used in weather, climate and ocean simulation systems [e.g. 8, 9, 10, respectively]. MPDATA is a solver for systems of advection equations of the following form: that describe evolution of a scalar field ψ transported by the fluid flow with velocity v. Quoting Numerical Recipes once more, development of methods to numerically solve such problems "is an art as much as a science" [1,Sec. 20.1], and MP-DATA is an example of the state-of-the art in this field. MP-DATA is designed to accurately solve equation (1) in an arbitrary number of dimensions assuring positive-definiteness of scalar field ψ and incurring small numerical diffusion. All relevant MPDATA formulae are given in the text but are presented without derivation or detailed discussion. For a recent review of MPDATA-based techniques see Smolarkiewicz [11,and references therein].
In this paper we introduce and discuss object-oriented implementations of an MPDATA-based two-dimensional (2D) advection equation solver written in C++11 (ISO/IEC 14882:2011), Python [13] and Fortran 2008 (ISO/IEC 1539-1:2010). In the following section we introduce the three implementations briefly describing the algorithm itself and discussing where and how the OOP techniques may be applied in its implementation. The syntax and nomenclature of OOP techniques are used without introduction, for an overview of OOP in context of C++, Python and Fortran, consult for example [15,Part II], [16,Chapter 5] and [17,Chapter 11], respectively. The third section of this paper covers performance evaluation of the three implementations. The fourth section covers discussion of the tradeoffs of the programming language choice. The fifth section closes the article with a brief summary.
Throughout the paper we present the three implementations by discussing source code listings which cover the entire program code. Subsections 2.1-2.6 describe all three implementations, while subsequent sections 2.7-2.12 cover discussion of C++ code only. The relevant parts of Python and Fortran codes do not differ significantly, and for readability reasons are presented in Appendix P and Appendix F, respectively.
The entire code is licensed under the terms of the GNU General Public License license version 3 [18].
All listings include line numbers printed to the left of the source code, with separate numbering for C++ (listings prefixed with C, black frame), listing C.0 (C++) 1 // code licensed under the terms of GNU GPL v3 2 // copyright holder: University of Warsaw which provide a convenient way of switching to different precision.
All codes are structured in a way allowing compilation of the code in exactly the same order as presented in the text within one source file, hence every Fortran listing contains definition of a separate module.

Array containers
Solution of equation (1) using MPDATA implies discretisation onto a grid of the ψ and the Courant number C = v · ∆t ∆x fields, where ∆t is the solver timestep and ∆x is the grid spacing.
Presented C++ implementation of MPDATA is built upon the Blitz++ library 1 . Blitz offers object-oriented representation of n-dimensional arrays, and array-valued mathematical expressions. In particular, it offers loop-free notation for array arithmetics that does not incur creation of intermediate temporary objects. Blitz++ is a header-only library 2 -to use it, it is enough to include the appropriate header file, and optionally expose the required classes to the present namespace: Here arr_t, rng_t and idx_t serve as alias identifiers and are introduced in order to shorten the code.
The power of Blitz++ comes from the ability to express array expressions as objects. In particular, it is possible to define a function that returns an array expression; i.e. not the resultant array, but an object representing a "recipe" defining the operations to be performed on the arguments. As a consequence, the return types of such functions become unintelligible. Luckily, the auto return type declaration from the C++11 standard allows to simplify the code significantly, even more if used through the following preprocessor macro: 1 Blitz++ is a C++ class library for scientific computing which uses the expression templates technique to achieve high performance, see http://sf.net/projects/blitz/ 2 Blitz++ requires linking with libblitz if debugging mode is used listing C.3 (C++) 8 #define return_macro(expr) \ 9 -> decltype(safeToReturn(expr)) \ 10

{ return safeToReturn(expr); }
The call to blitz::safeToReturn() function is included in order to ensure that all arrays involved in the expression being returned continue to exist in the caller scope. For example, definition of a function returning its array-valued argument doubled, reads: auto f(arr_t x) return_macro(2*x). This is the only preprocessor macro defined herein.
For the Python implementation of MPDATA the NumPy 3 package is used. In order to make the code compatible with both the standard CPython as well as the alternative PyPy implementation of Python [19], the Python code includes the following sequence of import statements: First, the PyPy's built-in NumPy implementation named numpypy is imported if applicable (i.e. if running PyPy), and the lazy evaluation mode is turned on through the set_invalidation(False) call. PyPy's lazy evaluation obtained with the help of a just-in-time compiler enables to achieve an analogous to Blitz++ temporary-array-free handling of arrayvalued expressions (see discussion in section 3). Second, to match the settings of C++ and Fortran compilers used herein, the NumPy package is instructed to ignore any floating-point errors, if such an option is available in the interpreter 4 . The above lines conclude all code modifications that needed to be added in order to run the code with PyPy.
Among the three considered languages only Fortran is equipped with built-in array handling facilities of practical use in high-performance computing. Therefore, there is no need for using an external package as with C++ and Python. Fortran array-handling features are not object-oriented, though.

Containers for sequences of arrays
As discussed above, discretisation in space of the scalar field ψ(x, y) into its ψ [i, j] grid representation requires floating-point array containers. In turn, discretisation in time requires a container class for storing sequences of such arrays, i.e. {ψ [n] , ψ [n+1] }. Similarly the components of the vector field C are in fact a {C [x] , C [y] } array sequence.
Using an additional array dimension to represent the sequence elements is not considered for two reasons. First, the 3  Consequently the last element of any such sequence may be accessed at index -1, the last but one at -2, and so on.
In the Python implementation the built-in tuple type is used to store sequences of NumPy arrays. Employment of negative indices for handling from-the-end addressing of elements is a built-in feature of all sequence containers in Python.
Fortran does not feature any built-in sequence container capable of storing arrays, hence a custom arrvec_t type is introduced:  The arr_t type is defined solely for the purpose of overcoming the limitation of lack of an array-of-arrays construct, and its only member field is a two-dimensional array. An array of arr_t is used hereinafter as a container for sequences of arrays.
The arrptr_t type is defined solely for the purpose of overcoming Fortran's limitation of not supporting allocatables of pointers. arrptr_t's single member field is a pointer to an instance of arr_t. Creating an allocatable of arrptr_t, instead of a multi-element pointer of arr_t, ensures automatic memory deallocation.
Type arrptr_t is used to implement the from-the-end addressing of elements in arrvec_t. The array data is stored in the arrs member field (of type arr_t). The at member field (of type arrptr_t) stores pointers to the elements of arrs. at has double the length of arrs and is initialised in a cyclic manner so that the -1 element of at points to the last element of arrs, and so on. Assuming psi is an instance of arrptr_t, the (i,j) element of the n-th array in psi may be accessed with psi%at( n )%p%a( i, j ).
The ctor(n) method initialises the container for a given number of elements n. The init(n,i,j) method initialises the n-th element of the container with a newly allocated 2D array spanning indices i(1):i(2), and j(1):j(2) in the first, and last dimensions respectively 6 . The so-called Arakawa-C staggered grid [20] depicted in Figure 1 is a natural choice for MPDATA. As a consequence, the discretised representations of the ψ scalar field, and each component of the C = v · ∆t ∆x vector field in eq. (1) are defined over different grid point locations. In mathematical notation this can be indicated by usage of fractional indices, e.g. C [x] [i− 1 /2, j] ,

Staggered grid
[i, j+ 1 /2] to depict the grid values of the C vector components surrounding ψ [i, j] . However, fractional indexing does not have a built-in counterpart in any of the employed programming languages. A desired syntax would translate i − 1 /2 to i − 1 and i + 1 /2 to i. OOP offers a convenient way to implement such notation by overloading the + and -operators for objects representing array indices.
In the C++ implementation first a global instance h of an empty structure hlf_t is defined, and then the plus and minus operators for hlf_t and rng_t are overloaded: This way, the arrays representing vector field components can be indexed using (i+h,j), (i-h,j) etc. where h represents the half.
In NumPy in order to prevent copying of array data during slicing one needs to operate on the so-called array views. Array views are obtained when indexing the arrays with objects of the Python's built-it slice type (or tuples of such objects in case of multi-dimensional arrays). Python forbids overloading of operators of built-in types such as slices, and does not define addition/subtraction operators for slice and int pairs. Consequently, a custom logic has to be defined not only for fractional indexing, but also for shifting the slices by integer intervals (i ± 1). It is implemented here by declaring a shift class with the adequate operator overloads: In Fortran fractional array indexing is obtained through definition and instantiation of an object representing the half, and having appropriate operator overloads: 7 One could argue that not using an own implementation of a slicerepresenting class in NumPy is a design flaw -being able to modify behaviour of a hypothetical numpy.slice class through inheritance would allow to implement the same behaviour as obtained in listing P.3 without the need to represent the unity as a separate object

Halo regions
The MPDATA formulae defining (discussed in the following sections) feature terms such as . One way of assuring validity of these formulae on the edges of the domain (e.g. for i=0) is to introduce the so-called halo region surrounding the domain. The method of populating the halo region with data depends on the boundary condition type. Employment of the halo-region logic implies repeated usage of array range extensions in the code such as i i ± halo. An ext() function is defined in all three implementation, in order to simplify coding of array range extensions:  Consequently, a range depicted by i ± 1/2 may be expressed in the code as ext(i, h). In all three implementations the ext() function accept the second argument to be an integer or a "half" (cf. section 2.3).

Array index permutations
Hereinafter, the π d a,b symbol is used to denote a cyclic permutation of an order d of a set {a, b}. It is used to generalise the MPDATA formulae into multiple dimensions using the following notation: Blitz++ ships with the RectDomain class (aliased here as idx_t) for specifying array ranges in multiple dimensions. The π permutation is implemented in C++ as a function pi() returning an instance of idx_t. In order to ensure compile-time evaluation, the permutation order is passed via the template parameter d (note the different order of i and j arguments in the two template specialisations): NumPy uses tuples of slices for addressing multidimensional array with a single object. Therefore, the following definition of function pi() suffices to represent π: listing P.6 (Python) 38 def pi(d, * idx): In the Fortran implementation pi() returns a pointer to the array elements specified by i and j interpreted as (i,j) or (j,i) depending on the value of the argument d. In addition to pi(), a helper span() function returning the length of one of the vectors passed as argument is defined: The span() function is used to shorten the declarations of arrays to be returned from functions in the Fortran implementation (see listings F.11 and F. 17-F.20).
It is worth noting here that the C++ implementation of pi() is branchless thanks to employment of template specialisation. With Fortran one needs to rely on compiler optimisations to eliminate the conditional expression within the pi() that depends on value of d which is always known at compile time.

Prototype solver
The tasks to be handled by a prototype advection equation solver proposed herein are: (i) storing arrays representing the ψ and C fields and any required housekeeping data, (ii) allocating and deallocating the required memory, (iii) providing access to the solver state, (iv) performing the integration by invoking the advectionoperator and boundary-condition handling routines.
In the following C++ definition of the solver structure, task (i) is represented with the definition of the structure member fields; task (ii) is split between the solver's constructor and the destructors of arrvec_t; task (iii) is handled by the accessor methods; task (iv) is handled within the solve method: The solver structure is an abstract definition (containing a pure virtual method) requiring its descendants to implement at least the advop() method which is expected to fill psi[n+1] with an updated (advected) values of psi[n]. The two template parameters bcx_t and bcy_t allow the solver to operate with any kind of boundary condition structures that fulfil the requirements implied by the calls to the methods of bcx and bcy, respectively.
The donor-cell and MPDATA schemes both require only the previous state of an advected field in order to advance the solution. Consequently, memory for two time levels (ψ [n] and ψ [n+1] ) is allocated in the constructor. The sizes of the arrays representing the two time levels of ψ are defined by the domain size (nx × ny) plus the halo region. The size of the halo region is an argument of the constructor. The cycle() method is used to swap the time levels without copying any data.
The arrays representing the C [x] and C [y] components of C, require (nx+1) × ny and nx × (ny+1) elements, respectively (being laid out on the Arakawa-C staggered grid).
Python definition of the solver class follows closely the C++ structure definition: The key difference stems from the fact that, unlike Blitz++, NumPy does not allow an array to have arbitrary index basein NumPy the first element is always addressed with 0. Consequently, while in C++ (and Fortran) the computational domain is chosen to start at (i=0, j=0) and hence a part of the halo region to have negative indices, in Python the halo region starts at (0,0) 8 . However, since the whole halo logic is hidden within the solver, such details are not exposed to the user. The bcx and bcy boundary-condition specifications are passed to the solver through constructor-like __init__() method as opposed to template parameters in C++.
The above C++ and Python prototype solvers in principle allow to operate with any boundary condition objects that implement methods called from within the solver. This requirement is checked at compile-time in the case of C++, and at run-time in the case of Python. In order to obtain an analogous behaviour with Fortran, it is required to define, prior to definition of a solver type, an abstract type with deferred procedures having abstract interfaces [sic!, see Table 2.1 in 21, for a summary of approximate correspondence of OOP nomenclature between Fortran and C++]: Having defined the abstract type for boundary-condition objects, a definition of a solver class following closely the C++ and Python counterparts may be provided:

Periodic boundaries (C++)
From this point, only C++ implementation is explained in the main text. The Python and Fortran implementations are included in appendices P and F.
The solver definition described in section 2.6 requires a given boundary condition object to implement a fill_halos() method. An implementation of periodic boundary conditions in C++ is provided in the following listing: a(pi<d>(left_halo, j)) = a(pi<d>(rght_edge, j)); 127 a(pi<d>(rght_halo, j)) = a(pi<d>(left_edge, j)); 128 } 129 }; As hinted by the member field names, the fill_halos() methods fill the left/right halo regions with data from the right/left edges of the domain. Thanks to employment of the function pi() described in section 2.5 the same code may be applied in any dimension (here being a template parameter).
Listings P.8 and F.8 contain the Python and Fortran counterparts to listing C.9.

Donor-cell formulae (C++)
MPDATA is an iterative algorithm in which each iteration takes the form of the so-called donor-cell formula (which itself is a first-order advection scheme).
MPDATA and donor-cell are explicit forward-in-time algorithms -they allow to predict ψ [n+1] as a function of ψ [n] where n and n + 1 denote two adjacent time levels. The donor-cell scheme may be written as [eq. 2 in 7]: where N is the number of dimensions, and F is the so-called flux function [7, eq. 3]: The flux function takes the following form in C++: Listings P.9-P11 and F.9-F.13 contain the Python and Fortran counterparts to listings C.12-C.15.

Donor-cell solver (C++)
As mentioned in the previous section, the donor-cell formula constitutes an advection scheme, hence we may use it to create a solver_donorcell implementation of the abstract solver class: The above definition is given as an example only. In the following sections an MPDATA solver of the same structure is defined.
Listings P.12 and F.14 contain the Python and Fortran counterparts to listing C.16.

MPDATA formulae (C++)
MPDATA introduces corrective steps to the algorithm defined by equation 2 and 3. Each corrective step is a donorcell step (eq. 2) with the Courant number fields corresponding to the MPDATA antidiffusive velocities of the following form [eqs 13, 14 in 7]: where ψ and C represent values from the previous iteration and where: For positive-definite ψ, the A and B terms take the following form 9 : If the denominator in equations 6 or 7 equals zero for a given i and j, the corresponding A [i, j] and B [i, j] are set to zero what may be conveniently represented with the where construct (available in all three considered languages):

MPDATA solver (C++)
An MPDATA solver may be now constructed by inheriting from solver class with the following definition in C++: The array of sequences of temporary arrays tmp allocated in the constructor is used to store the antidiffusive velocities from the present and optionally previous timestep (if using more than two iterations).
The advop() method controlls the MPDATA iterations within one timestep. The first (step = 0) iteration of MPDATA is an unmodified donor-cell step (compare listing C.15). Subsequent iterations begin with calculation of the antidiffusive Courant fields using formula 4. In order to calculate values spanning an (i− 1 ⁄2 ... i+ 1 ⁄2) range using a formula for C [i+ 1 /2,...] only, the formula is evaluated using extended index ranges im and jm. In the second (step=1) iteration the uncorrected Courant field (C_unco) points to the original C field, and the antidiffusive Courant field is written into C_corr which points to tmp [

Usage example (C++)
The following listing provides an example of how the MP-DATA solver defined in section 2.11 may be used together with the cyclic boundary conditions defined in section 2.7. In the example a Gaussian signal is advected in a 2D domain defined over a grid of 24×24 cells. The program first plots the initial condition, then performs the integration for 75 timesteps with three different settings of the number of iterations used in MPDATA. The velocity field is constant in time and space (although it is not assumed in the presented implementations). The signal shape at the end of each simulation is plotted as well. Plotting is done with the help of the gnuplot-iostream library 10 .
The resultant plot is presented herein as Figure 2. The top panel depicts the initial condition. The three other panels show a snapshot of the field after 75 timesteps. The donor-cell solution is characterised by strongest numerical diffusion resulting in significant drop in the signal amplitude. The signals advected using MPDATA show smaller numerical diffusion with the solution obtained with more iterations preserving the signal altitude more accurately. In all of the simulations the signal maintains its positive definiteness. The domain periodicity is apparent in the plots as the maximum of the signal after 75 timesteps is located near the domain walls.

Performance evaluation
The three introduced implementations of MPDATA were tested with the following set-ups employing free and opensource tools:

C++:
• GCC g++ 4.8.0 11  The performance tests were run on a Debian and an Ubuntu GNU/Linux systems with the above-listed software obtained via binary packages from the distributions' package repositories (most recent package versions at the time of writing). The tests were performed on two 64-bit machines equipped with an AMD Phenom ™ II X6 1055T (800 MHz) and an Intel ® Core ™ i5-2467M (1.6 GHz) processors.
For both C++ and Fortran the GCC compilers were invoked with the -Ofast and the -march=native options. The Clang compiler was invoked with the -O3, the -mllvm -vectorize, the -ffast-math and the -march=native options. The CPython interpreter was invoked with the -OO option.
In addition to the standard Python implementation CPython, the Python code was tested with PyPy. PyPy is an alternative implementation of Python featuring a just-in-time compiler. PyPy includes an experimental partial reimplementation of NumPy that compiles NumPy expressions into native assembler. Thanks to employment of lazy evaluation of array expressions (cf. Sect. 2.1) PyPy allows to eliminate the use of temporary matrices for storing intermediate results, and to perform multiple operations on the arrays within a single array index traversal 12 . Consequently, PyPy allows to overcome the same performance-limiting factors as those addressed by Blitz++, although the underlying mechanisms are different. In contrast to other solutions for improving performance of NumPy-based codes such as Cython 13 , numexpr 14 or Numba 15 , PyPy does not require any modifications to the code. Thus, PyPy may serve as a drop-in replacement for CPython ready to be used with previously-developed codes.
The same set of tests was run with all four set-ups. Each test set consisted of 16 program runs. The test programs are analogous to the example code presented in section 2.12. The 11 GNU Compiler Collection packaged in the Debian's gcc-snapshot_20130222-1 12 Lazy evaluation available in PyPy 1.9 has been temporarily removed from PyPy during a refactoring of the code. It'll be reinstantiated in the codebase as soon as possible, but past PyPy 2.0 release 13   tests were run with different grid sizes ranging from 64×64 to 2048×2048. The Gaussian impulse was advected for nt = 2 24 /(nx · ny) timesteps (2 24 chosen arbitrarily), in order to assure comparable timing accuracy for all grid sizes. Three MP-DATA iterations were used (i.e. two corrective steps). The initial condition was loaded from a text file, and the final values were compared at the end of the test with values loaded from another text file assuring the same results were obtained with all four set-ups. The tests were run multiple times; program startup, data loading, and output verification times were subtracted from the reported values (see caption of Figure 4 for details). • Fortran gives shortest execution times for any domain size; • C++ execution times are less than twice those of Fortran for grids larger than 256×256; • CPython requires from around 4 to almost 10 times more CPU time than Fortran depending on the grid size; • PyPy execution times are in most cases closer to C++ than to CPython. The support for OOP features in gfortran, the NumPy support in PyPy, and the relevant optimisation mechanisms in GCC are still in active development and hence the performance with some of the set-ups may likely change with newer versions of these packages.
It is worth mentioning, that even though the three implementations are equally structured, the three considered languages have some inherent differences influencing the execution times. Notably, while Fortran and Blitz++ offer runtime array-bounds and array-shape checks as options not intended for use in production binaries, NumPy performs them always. Additionally, the C++ and Fortran set-ups may, in principle, benefit from GCC's auto-vectorisation features which do not have yet counterparts in CPython or PyPy. Finally, Fortran uses different ordering for storing array elements in memory, but since all tests were carried out using square grids, this should not have had any impact on the performance 17 .
The authors do expect some performance gain could be obtained by introducing into the codes some "manual" optimisations -code rearrangements aimed solely at the purpose of increasing performance. These were avoided intentionally as they degrade code readability, should in principle be handled by the compilers, and are generally advised to be avoided [e.g. 22, section 3.12].

Discussion on the tradeoffs of language choice
One of the aims of this paper is to show the applicability of OOP features of the three programming languages (or language-library pairs) for scientific computing. The main focus is to represent what can be referred to as blackboard abstractions [21] within the code. Presented benchmark tests, although quite simplistic, together with the experience gained from the development of codes in three different languages provide a basis for discussion on the tradeoffs of programming language choice. The discussion concerns in principle the development of finite-difference solvers for partial differential equations, but is likely applicable to the scientific software in general. A partly objective and partly subjective summary of pros and cons of C++, Python and Fortran is presented in the four following subsections.

OOP for blackboard abstractions
It was shown in section 2 that C++11/Blitz++, Python/NumPy and Fortran 2008 provide comparable functionalities in terms of matching the blackboard abstractions within the program code. Taking into account solely the part of code representing particular formulae (e.g. listings C.21, P.17, F.20 and equation 4) all three languages allow to match (or surpass) L A T E X in its brevity of formula translation syntax. All three languages were shown to be capable of providing mechanisms to compactly represent such abstractions as: • loop-free array arithmetics; • definitions of functions returning array-valued expressions; • permutations of array indices allowing dimensionindependent definitions of functions (see e.g. listings C.12 and C.13, P.10 and P.11, F.11 and F.12); • fractional indexing of arrays corresponding to employment of a staggered grid. Three issues specific to Fortran that resulted in employment of a more repetitive or cumbersome syntax than in C++ or Python were observed: • Fortran does not feature a mechanism allowing to reuse a single piece of code (algorithm) with different data types (compare e.g. listings C.6, P.5 and F.4) such as templates in C++ and the so-called duck typing in Python; • Fortran does not allow function calls to appear on the left hand side of assignment (see e.g. how the ptr pointers were used as a workaround in the cyclic_fill_halos method in listing F.8); • Fortran lacks support for arrays of arrays (cf. Sect. 2.2). Interestingly, the limitation in extendability via inheritance was found to exist partially in NumPy as well (see footnote 7). The lack of a counterpart in Fortran to the C++ template mechanism was identified in [23] as one of the key deficiencies of Fortran when compared with C++ in context of applicability to objectoriented scientific programming.

Performance
The timing and memory usage statistics presented in figures 3-5 reveal that no single language/library/compiler set-up corresponded to both shortest execution time and smallest memory footprint.
One may consider performance measures addressing not only the program efficiency but also the factors influencing the development and maintenance time/cost [of particular importance in scientific computing, 24]. Taking into account such measures as code length or coding time, the Python environment gains significantly. Presented Python code is shorter than the C++ and Fortran counterparts, and is simpler in terms of syntax and usage (see discussion below).
Employment of the PyPy drop-in replacement for the standard Python implementation brings Python's performance significantly closer to those of C++ and Fortran, in some cases making it the least memory consuming set-up. Python has already been the language of choice for scientific software projects having code clarity or ease of use as the first requirement [see e.g. 25]. PyPy's capability to improve performance of unmodified Python code may make Python a favourable choice even if high performance is important, especially if a combined measure of performance and development cost is to be considered.

Ease of use and abuse
Using the number of lines of code or the number of distinct language keywords needed to implement the MPDATA-based solver presented in section 2 as measures of syntax brevity, Python clearly surpasses its rivals. Python was developed with emphasis on code readability and object-orientation. Arguably, taking it to the extreme -Python uses line indentation to define blocks of code and treats even single integers as objects. As a consequence Python is easy to learn and easy to teach. It is also much harder to abuse Python than C++ or Fortran (for instance with goto statements, employment of the preprocessor, or the implicit typing in Fortran).
Python implementations do not expose to the user the compilation or linking processes. As a result, Python-written software is easier to deploy and share, especially if multiple architectures and operating systems are targeted. However, there exist tools such as CMake 18 that allow to efficiently automate building, testing and packaging of C++ and Fortran programs.
Python is definitely easiest to debug among the three languages. Great debugging tools for C++ do exist, however the debugging and development is often hindered by indecipherable compiler messages flooded with lengthy type names stemming from employment of templates. Support for the OOP features of Fortran among free and open source compilers, debuggers and other programming aids remains immature.
With both Fortran and Python, the memory footprint caused by employment of temporary objects in array arithmetics is dependant on compiler choice or the level of optimisations. In contrast, Blitz++ ensures temporary-array-free computations by design [26] avoiding unintentional performance loss.

Added values
The size of the programmers' community of a given language influences the availability of trained personnel, reusable software components and information resources. It also affects the maturity and quality of compilers and tools. Fortran is a domain-specific language while Python and C++ are generalpurpose languages with disproportionately larger users' communities. The OOP features of Fortran have not gained wide popularity among users [27] 19 . Fortran is no longer routinely taught at the universities [28], in contrast to C++ and Python. An example of decreasing popularity of Fortran in academia is the discontinuation of Fortran printed editions of the "Numerical Recipes" series of Press et al.
Blitz++ is one of several packages that offer highperformance object-oriented array manipulation functionality with C++ (and is not necessarily optimal for every purpose [29]). In contrast, the NumPy package became a de facto standard solution for Python. Consequently, numerous Python libraries adopted NumPy but there are apparently very few C++ libraries offering Blitz++ support out of the box (the gnuplotiostream used in listing C.24 being a much-appreciated counterexample). However, Blitz++ allows to interface with virtually any library (including Fortran libraries), by resorting to referencing the underlying memory with raw pointers.
The availability and quality of libraries that offer objectoriented interfaces differs among the three considered languages. The built-in standard libraries of Python and C++ are richer than those of Fortran and offer versatile data types, collections of algorithms and facilities for interaction with host operating system. In the authors' experience, the small popularity of OOP techniques among Fortran users is reflected in the library designs (including the Fortran's built-in library routines). What makes correct use of external libraries more difficult with Fortran is the lack of standard exception handling mechanism, a feature long and much requested by the numerical community [30,Foreword].
Finally, the three languages differ as well with regard to availability of mechanisms for leveraging shared-memory parallelisation (e.g. with multi-core processors). GCC supports OpenMP with Fortran and C++. The CPython and PyPy implementations of Python do not offer any built-in solution for multi-threading.

Summary and outlook
Three implementations of a prototype solver for the advection equation were introduced. The solvers are based on MP-DATA -an algorithm of particular applicability in geophysical fluid dynamics [11]. All implementations follow the same object-oriented structure but are implemented in three different languages: • C++ with Blitz++; • Python with NumPy; • Fortran. Presented programs were developed making use of such recent developments as support for C++11 and Fortran 2008 in GCC, and the NumPy support in the PyPy implementation of Python. The fact that all considered standards are open and the employed tools implementing them are free and open-source is certainly an advantage [31].
The key conclusion is that all considered language/library/compiler set-ups offer possibilities for using OOP to compactly represent the mathematical abstractions within the program code. This creates the potential to improve code readability and brevity, • contributing to its auditability, indispensable for credible and reproducible research in computational science [32,33,34]; and • helping to keep the programs maintainable and avoiding accumulation of the code debt 20 that besets scientific software in such domains as climate modelling [36]. The performance evaluation revealed that: • the Fortran set-up offered shortest execution times, • it took the C++ set-up less than twice longer to compute than Fortran, • C++ and PyPy set-ups offered significantly smaller memory consumption than Fortran and CPython for larger domains, • the PyPy set-up was roughly twice slower than C++ and up to twice faster than CPython. The three equally-structured implementations required ca. 200, 300, and 500 lines of code in Python, C++ and Fortran, respectively.
In addition to the source code presented within the text, a set of tests and build-/test-automation scripts allowing to reproduce the analysis and plots presented in section 3 are all available in the CPC Program Library and at the project repository 21 , and are released under the GNU GPL license [18]. The authors encourage to use the presented codes for teaching and benchmarking purposes.
The OOP design enhances the possibilities to reuse and extend the presented code. Development is underway of an object-oriented C++ library featuring concepts presented herein, supporting integration in one to three dimensions, handling systems of equations with source terms, providing miscellaneous options of MPDATA and several parallel processing approaches.