Regrid
The mjolnir library/folder contains tools to transform from the icosahedral grid (not convenient for plotting) to a latitude-longitude grid (easy for plotting). The regrid command can be invoked on the command line, similar to invoking mjolnir:
regrid -i <initial_file_index> -l <last_file_index> <simulation_folder>
So, for example
regrid -i 10 -l 20 earth_hs
will regrid the data in the folder earth_hs from file index 10 to 20.
The current regrid process requires PyCUDA and a GPU. There is some old code in the file hamarr_suppl.py that does the regridding using a CPU and various other techniques, but you are on your own if you decide to use it.
For some types of plots, it makes sense use pressure as the vertical coordinate. This requires doing some interpolation from altitude/height to pressure. Currently, we do this by simple linear interpolation (old versions of the mjolnir/regrid code used PCHIP interpolation, which is almost certainly overkill). The first tricky part of this is: how do you decide what the pressure coordinates should be, given that the pressure varies throughout the simulation? The other tricky part: if you are doing time-averages, the pressure grid across all points in time needs to be the same, unless you are doing further interpolation.
I've done my best to supply a flexible solution to these issues. There is a function in the mjolnir libraries that reads in the pressure data from a file or group of files, and averages the pressure in each altitude layer, then averages that in time (if more that one file is requested), to produce a pgrid text file containing a list of pressure that will then be used as the new vertical coordinate. The naming convention is pgrid_x_y_z.txt, where x is <initial_file_index>, y is <last_file_index>, and z is the stride between x and y.
The pgrid function can be invoked by itself from the command line, similar to regrid:
pgrid -i <initial_file_index> -l <last_file_index> <simulation_folder> -stride <stride_from_x_to_y>
(-stride defaults to 1, so that every file from <initial_file_index> to <last_file_index> is included).
When running regrid or mjolnir, these functions will autogenerate any pgrid files needed, if the -pgrid option is auto (which is also the default value, if it is omitted). The auto behavior is to pass <initial_file_index> and <last_file_index> from regrid or mjolnir to pgrid, with the stride set to 1. One can also specify an existing filename for the -pgrid option, if a specific vertical grid is desired.
The regrid process then emits two output files for every output time. One will be a latitude-longitude-height grid inside <simulation_folder>/, with a naming convention like regrid_height_<planet>_<index>.h5. The other will be contained in <simulation_folder>/<pgrid_file>/' with a name like regrid__.h5`. This means that you can have many different regrids in separate folders, all using a different pressure grid.
In running mjolnir to produce plots, with -vc pressure (default) and -pgrid auto (default), the code will look for a pgrid file and folder that matching the default behavior described above (with x = <initial_file_index>, y = <last_file_index>, and z = 1). If it does not exist, it runs pgrid and regrid to generate those files. As with regrid, one can also set -pgrid to a specific file, if that file and its corresponding folder with regridded data exist.
For mapping the height grid onto a pressure grid, we use simple linear interpolation/extrapolation (see above). If the pressure does not decrease monotonically with height, the code will raise an error.
Mapping the icosahedral grid onto a latitude-longitude grid is somewhat more complicated. To do this, we utilize the Möller-Trumbore algorithm. The basic idea is to calculate the intersection of a ray (represented by the destination latitude-longitude point) and a triangle (the three nearest points of the icosahedral grid). The position of the ray can is written as a function of the locations of the three vertices of the triangle and three weights. We can then construct a matrix of these weights, which then operates on a matrix of the field values on the icosahedral grid to produce a matrix of the field values on the desired latitute-longitude grid.
In using the Möller-Trumbore algorithm, we are treating each triangle (formed by three nearby points) as if it were flat (not part of a sphere), but the approximation is typically a good one. In any case, it produces results that are much smoother than a nearest neighbors interpolation and is free from the boundary issues one would have at the poles and across the 360 degree longitude line, if the rectilinear interpolation is used.
The horizontal resolution of the destination grid is taken to be approximately the same as the original grid. Specifically, the angular size of grid boxes is given (in degrees) by:
The parameter g_level controls the horizontal resolution of the original icosahedral grid.
PyCUDA is used to parallelize both the vertical and horizontal regridding as much as possible.