Posted on 23 July 2014
[This is a longish post about the process of developing pyrrtm, a python wrapper for RRTM. Please see the project page for more practical details about using it!]
The RRTM radiation code is used a lot in atmospheric science, and is pretty powerful. Unfortunately, it is stuck with a rather nasty ASCII interface (which I can’t really imagine anyone using!) so when I came to use it to do some work, I wrote a python wrapper. I actually went through various iterations for how this wrapper should work, and initially I considered writing a script that could construct the ASCII input required by RRTM. I quickly ruled this out as it would be pretty tedious, and also introduces a large overhead. We ended up doing ~10,000–100,000 radiative transfer calculations for each result in the paper, so overhead was important for us.
After ruling out using the ASCII interface, the first thing I tried was using the excellent f2py to compile the FORTRAN sources for RRTM into a python module. This proved successful and surprisingly easy (once I’d worked out which RRTM variables needed setting!), as well as having no overhead at all — numpy arrays are fed straight into the RRTM code. This was the method I used to produce the results for the paper. After completing that project, I intended to release the code, but I hit an obstacle. For some reason, I had difficulties compiling the code, with segfaults on new versions of gfortran, and unpredictable behaviour (i.e. sometimes I got the wrong answer with no error). I put this down to the extremely old-school nature of the RRTM library and sensitivity to compiler options. With f2py, it was too difficult to have control over those things. Sadly, this unpredictable behaviour meant that releasing the code was a bit of a waste of time. :-(
After some time doing other projects, I eventually got around to attempting a solution to these issues. First, I added a little NetCDF wrapper (written in C) to the RRTM code. I simplified matters by only exposing the parts of RRTM we needed for our project (and similar projects). As the wrapper is compiled into an application with RRTM, I had something that always worked as a standalone application, and with relatively straightforward compiler flags (as it did not need to be a shared object and link with python). I also fixed a couple of bugs caused by odd gfortran behaviour hitting up against some rather dodgy FORTRAN (issues with DOUBLE PRECISION declarations clashing with -fdefault-real-8 flag — really difficult to fix!). Next, I redid the python module interface using cython in place of f2py to give better control over compilation. This was very simple, and on reflection I now think that cython is better than f2py for exposing small sets of FORTRAN functions to python.
I now have two ways to interface with RRTM (the cython module and NetCDF), both of which are better than the original ASCII interface. I still occasionally have issues getting the cython module approach to work, but the NetCDF approach seems rock solid! If speed is critical, the cython module wins (no overhead – up to 10 times speed up on my machine), but for other purposes the NetCDF method is the best. Both methods are functionally identical, and just vary in performance really.
This section is now out of date! The cython wrapper now raises exceptions :-)
After battling to produce not one but two python interfaces, the easy part was writing a bit of python to call the methods! The python package, pyrrtm, can toggle between the two methods, letting the user choose. Here are the pros and cons of each method once wrapped by pyrrtm (at time of writing!):
| Pros | Cons | |
|---|---|---|
| NetCDF |
|
|
| cython |
|
|
As you can see, the cython interface has the unfortunate disadvantage of killing python when it fails rather than gracefully raising an exception. This means that the NetCDF interface is always default in pyrrtm, but the cython interface can easily be switched on whenever needed (e.g. you can test using the NetCDF interface, but use cython for a big job).
Overall I’m pretty happy with how it turned out, and please contact me if you plan on using it or want to contribute to it!