Fortran & Python Integration

The two-language problem

The ecosystem lesson's history lesson repeated itself here: scientists prototype in Python (readable, dynamic, rich libraries), then hit a loop that takes hours, rewrite it in Fortran, and lose weeks gluing the two. The mature pattern flips the balance: Python orchestrates, Fortran computes. Python owns loops that orchestrate — file walking, experiment control, plotting — seconds-scale code where readability wins. Fortran owns the inner numerical kernels — the per-time-step physics — where a 50× speedup is routine. The integration layer is the only piece worth engineering.

Level 1: subprocess and files

The simplest bridge: compile a standalone Fortran executable, launch it from Python with subprocess.run, exchange data through files. Zero libraries, zero compilation toolchain on the Python side, and trivially debuggable — each side is an ordinary program. The cost is serialization and process startup, so this is prototyping glue, not production HPC glue. Everything the files lesson taught about formats is exactly what this level consumes.

import subprocess

# run the Fortran solver once; it reads config from stdin, writes CSV out
subprocess.run(["./solver", "case_07.nml"], check=True)

# the Python side now parses the produced file like any data source
with open("case_07.csv") as f:
    for line in f:
        ...

Fig. 1 — Files first, then the zero-copy extension, then the raw C ABI — each level trades convenience for speed.

Level 2: F2PY — the production bridge

F2PY — part of NumPy since its beginnings — reads your Fortran source, generates a C/CPython extension, and compiles it into a module Python imports directly. Fortran subroutines become Python functions; Fortran arrays become NumPy arrays pointing at the same memory — no copy, one of the few true zero-copy bridges in the scientific stack. The Fortran side needs nothing special except intent annotations, which F2PY uses to decide which arguments are outputs:

! solver.f90 — the compute kernel, ordinary Fortran
subroutine euler_step(n, dt, y, ynext)
  ! One time step of a simple explicit integrator.
  integer, intent(in) :: n
  real(8), intent(in) :: dt, y(n)
  real(8), intent(out) :: ynext(n)
  ynext = y + dt * sin(y)          ! whole-array math, vectorizes
end subroutine euler_step

# build the extension module (Windows: replace gfortran with your mingw setup)
python -m numpy.f2py -c solver.f90 -m solver

import numpy as np
import solver                     # our Fortran, now a Python module

y = np.zeros(10_000)
dt = 0.01
for step in range(1000):          # the ORCHESTRATION loop stays in Python
    y = solver.euler_step(y, dt)  # the COMPUTE call lands in Fortran
print(y.sum())

Two details make it work smoothly. F2PY wraps intent(out) array arguments by creating the result array and returning it (as shown), and it accepts intent(in, out) for in-place updates. Shape information flows both ways: pass a (n,) NumPy array and n can be inferred from its shape if the dummy is declared with an assumed size or explicit shape using another dummy. The modern f2py==2.x uses meson builds and handles most iso_c_binding interfaces automatically.

Level 3: ctypes and the C ABI

When the Fortran side exposes bind(c) symbols (interop lesson), any tool that calls C — ctypes, CFFI, ctypes from Rust, Java's JNA — can call Fortran directly. This is the "shared library as API" route: gfortran -shared -fPIC libsolver.f90 -o libsolver.so, then Python declares argument types explicitly and calls through the C function-pointer table. No build-time coupling to NumPy, at the price of manual type declarations and no automatic output-array creation:

import ctypes
import numpy as np

lib = ctypes.CDLL("./libsolver.so")
lib.scale_array.argtypes = [ctypes.c_int,
                            np.ctypeslib.ndpointer(dtype=np.float64,
                                                   ndim=1, flags="C_CONTIGUOUS"),
                            ctypes.c_double]
x = np.ones(5)
lib.scale_array(x.size, x, 3.0)   # Fortran multiplies x in place
print(x)                          # [3. 3. 3. 3. 3.]

The array must be C-contiguous (the default for newly created NumPy arrays) because the Fortran side sees a raw address plus a length — this is the rawest level and therefore the one to use when many languages must share one library.

Design guidelines that survive contact

• Keep the hot loop in Fortran. The Python loop in the F2PY example calls Fortran 1000 times per second — that is the correct shape. The wrong shape calls Python from a Fortran loop, or copies data in and out each call.
• Batch, don't stream. If your Python counterpart needs every intermediate value, restructure: compute the whole trajectory in Fortran, return it once. Transfer time is amortized by the kernel's arithmetic.
• Validate once. Check array shapes and finite-ness on the Python side before the first call; the Fortran side assumes the contract (its intent is checked at its own compile time, not at the boundary).
• Debug the Fortran side with its own tools (gfortran -fcheck=bounds) before suspecting the bridge — the error lesson's backtraces print the Fortran call chain, not Python frames.