In atmospheric sciences lots of Fortran is used, and in-fact lots of Fortran is still being written. One day I wanted to validate some model code I had written against an older Fortran implementation, and instead of writing the 10 lines of Fortran I would have needed to write I decided to try and wrap the code in Python instead. How hard could it be? SciPy itself is 13.5% Fortran code, surely this is a solved problem. I could not have been more wrong.

I had heard of F2PY before and thought that would be the best place to start. The code that I was trying to wrap consisted was F90 style with a large number of modules and files, and I quickly realized that the best way to wrap it would be to separate it into it’s own library and write a Fortran interface to it and wrap that. Things went really downhill when I realized that F2PY cannot handle Fortran derived types.

If you have not had the privilege of working with Fortran bnfore, derived types are similar to a C structure. Here is an example from the code I was wrapping,

TYPE LIDORT_Fixed_Chapman
      REAL(fpk), dimension (0:MAXLAYERS) :: TS_HEIGHT_GRID
      REAL(fpk), dimension (0:MAXLAYERS) :: TS_PRESSURE_GRID
      REAL(fpk), dimension (0:MAXLAYERS) :: TS_TEMPERATURE_GRID
      INTEGER,   dimension (MAXLAYERS)   :: TS_FINEGRID ( MAXLAYERS )
      REAL(fpk) :: TS_RFINDEX_PARAMETER
END TYPE LIDORT_Fixed_Chapman

Since the main input to the model was ~10 of these derived types, I had to rethink my strategy. There were a few options:

  • I found a package called f90wrap which can bind derived types, and seemed very promising. The downside is that it requires a runtime component rather than being a static binding generator. It was unclear if it would work on my M1 macbook as well so I ultimately decided against it.

  • Fortran 2003 introduced iso_c_binding which allows for enhanced C interoperability. Essentially the code could be modified to make the derived types interoperable with C structures. The code could then be called in Python through either cffi or cython or something similar. I didn’t go this route because my Fortran knowledge is fairly poor, and it would require modifying the original source code in some way.

  • The third and simplest method is to write a wrapper function that that takes in all of the input and output values, constructs the derived types, and calls the model. This wrapper can then be used in the F2PY generator, and this is how I ended up with this monstrosity

subroutine internal_lidort(TS_DO_FULLRAD_MODE, TS_DO_THERMAL_EMISSION, TS_DO_SURFACE_EMISSION, TS_DO_PLANE_PARALLEL, &
                  TS_DO_BRDF_SURFACE, TS_DO_UPWELLING, TS_DO_DNWELLING, TS_DO_TOA_CONTRIBS, TS_DO_SURFACE_LEAVING, &
                  TS_DO_SL_ISOTROPIC, TS_DO_WATER_LEAVING, TS_DO_FLUORESCENCE, TS_DO_TF_ITERATION, &
                  TS_DO_WLADJUSTED_OUTPUT, TS_DO_TOA_ILLUMINATION, TS_DO_BOA_ILLUMINATION, TS_DO_ALBTRN_MEDIA, &
                  TS_DO_PLANETARY_PROBLEM, TS_DO_MSSTS, &
                  TS_TAYLOR_ORDER, TS_NSTREAMS, TS_NLAYERS, TS_NFINELAYERS, TS_N_THERMAL_COEFFS, TS_LIDORT_ACCURACY, &
                  TS_ASYMTX_TOLERANCE, TS_TF_MAXITER, TS_TF_CRITERION, TS_TOA_ILLUMINATION, TS_BOA_ILLUMINATION, TS_FLUX_FACTOR, &
                  TS_HEIGHT_GRID, TS_PRESSURE_GRID, TS_TEMPERATURE_GRID, TS_FINEGRID, TS_RFINDEX_PARAMETER, &
                  TS_DELTAU_VERT_INPUT, TS_PHASMOMS_TOTAL_INPUT, TS_PHASFUNC_INPUT_UP, TS_PHASFUNC_INPUT_DN, TS_LAMBERTIAN_ALBEDO, &
                  TS_THERMAL_BB_INPUT, TS_SURFACE_BB_INPUT, TS_ATMOS_WAVELENGTH, &
                  TS_DO_DEBUG_WRITE, TS_DO_WRITE_INPUT, TS_INPUT_WRITE_FILENAME, TS_DO_WRITE_SCENARIO, TS_SCENARIO_WRITE_FILENAME, &
                  TS_DO_WRITE_FOURIER, TS_FOURIER_WRITE_FILENAME, TS_DO_WRITE_RESULTS, TS_RESULTS_WRITE_FILENAME, &
                  TS_DO_FOCORR, TS_DO_FOCORR_EXTERNAL, TS_DO_FOCORR_NADIR, TS_DO_FOCORR_OUTGOING, TS_DO_SSCORR_TRUNCATION, &
                  TS_DO_SSCORR_USEPHASFUNC, &
                  TS_DO_EXTERNAL_WLEAVE, TS_DO_DOUBLE_CONVTEST, TS_DO_SOLAR_SOURCES, TS_DO_REFRACTIVE_GEOMETRY, &
                  TS_DO_CHAPMAN_FUNCTION, TS_N_USER_LEVELS, &
                  TS_DO_RAYLEIGH_ONLY, &
                  TS_DO_ISOTROPIC_ONLY, TS_DO_NO_AZIMUTH, TS_DO_ALL_FOURIER, TS_DO_DELTAM_SCALING, TS_DO_SOLUTION_SAVING, &
                  TS_DO_BVP_TELESCOPING, &
                  TS_DO_USER_STREAMS, TS_DO_ADDITIONAL_MVOUT, TS_DO_MVOUT_ONLY, TS_DO_THERMAL_TRANSONLY, &
                  TS_DO_OBSERVATION_GEOMETRY, &
                  TS_DO_DOUBLET_GEOMETRY, &
                  TS_NMOMENTS_INPUT, TS_NBEAMS, TS_BEAM_SZAS, TS_N_USER_RELAZMS, TS_USER_RELAZMS, &
                  TS_N_USER_STREAMS, TS_USER_ANGLES_INPUT, &
                  TS_USER_LEVELS, TS_GEOMETRY_SPECHEIGHT, TS_N_USER_OBSGEOMS, TS_USER_OBSGEOMS_INPUT, &
                  TS_N_USER_DOUBLETS, TS_USER_DOUBLETS, &
                  TS_EARTH_RADIUS, TS_OMEGA_TOTAL_INPUT, TS_N_SURFACE_WFS, &
                  TS_LAYER_VARY_FLAG, TS_LAYER_VARY_NUMBER, TS_N_TOTALCOLUMN_WFS, TS_N_SLEAVE_WFS, &
                  TS_COLUMNWF_NAMES, TS_PROFILEWF_NAMES, TS_L_DELTAU_VERT_INPUT, TS_L_OMEGA_TOTAL_INPUT, &
                  TS_L_PHASMOMS_TOTAL_INPUT, TS_L_PHASFUNC_INPUT_UP, TS_L_PHASFUNC_INPUT_DN, &
                  TS_DO_COLUMN_LINEARIZATION, TS_DO_PROFILE_LINEARIZATION, TS_DO_ATMOS_LINEARIZATION, &
                  TS_DO_SURFACE_LINEARIZATION, TS_DO_LINEARIZATION, TS_DO_SIMULATION_ONLY, TS_DO_ATMOS_LBBF, &
                  TS_DO_SURFACE_LBBF, TS_DO_SLEAVE_WFS, NUM_REPEAT, &
                  TS_INTENSITY, TS_MEANI_DIFFUSE, TS_FLUX_DIFFUSE, TS_DNMEANI_DIRECT, TS_DNFLUX_DIRECT, &
                  TS_ALBMED_USER, TS_TRNMED_USER, TS_ALBMED_FLUXES, TS_TRNMED_FLUXES, TS_PLANETARY_TRANSTERM, &
                  TS_PLANETARY_SBTERM, TS_PATHGEOMS, TS_LOSTRANS, TS_LAYER_MSSTS, TS_SURF_MSSTS, &
                  TS_CONTRIBS, TS_COLUMNWF, TS_MEANI_DIFFUSE_COLWF, TS_FLUX_DIFFUSE_COLWF, TS_DNMEANI_DIRECT_COLWF, &
                  TS_DNFLUX_DIRECT_COLWF, &
                  TS_PROFILEWF, TS_MEANI_DIFFUSE_PROFWF, TS_FLUX_DIFFUSE_PROFWF, TS_DNMEANI_DIRECT_PROFWF, &
                  TS_DNFLUX_DIRECT_PROFWF, &
                  TS_ABBWFS_JACOBIANS, TS_ABBWFS_FLUXES, TS_ALBMED_USER_PROFWF, TS_TRNMED_USER_PROFWF, &
                  TS_ALBMED_FLUXES_PROFWF, &
                  TS_TRNMED_FLUXES_PROFWF, TS_TRANSBEAM_PROFWF, TS_ALBMED_USER_COLWF, TS_TRNMED_USER_COLWF, &
                  TS_ALBMED_FLUXES_COLWF, &
                  TS_TRNMED_FLUXES_COLWF, TS_TRANSBEAM_COLWF, TS_PLANETARY_TRANSTERM_PROFWF, TS_PLANETARY_SBTERM_PROFWF, &
                  TS_PLANETARY_TRANSTERM_COLWF, TS_PLANETARY_SBTERM_COLWF, &
                  TS_LC_LOSTRANS, &
                  TS_LC_LAYER_MSSTS, TS_LC_SURF_MSSTS, TS_LP_LOSTRANS, TS_LP_LAYER_MSSTS, TS_LP_SURF_MSSTS, &
                  TS_SURFACEWF, &
                  TS_MEANI_DIFFUSE_SURFWF, TS_FLUX_DIFFUSE_SURFWF, TS_SBBWFS_JACOBIANS, TS_SBBWFS_FLUXES, TS_LS_LAYER_MSSTS, &
                  TS_LS_SURF_MSSTS)

It did work, by the time I realized how ridiculous this was going to be I was too far lost in my madness. If I could go back and re-do it I would have learned how to use iso_c_binding instead and abandoned F2PY. While this solution works, it probably involves a lot of unnecessary memory copies when constructing the derived types.

What I have learned is that if you are trying to wrap older, F77 code, or code that does not heavily use derived types in it’s interface, you can use F2PY and be fine. If you want to wrap a more complicated Fortran model, then prepare to be in for a bit of a headache.