Fortran Is Not Dead

June 3, 2026 20 min read Software Development

fortranhpcscientific-computingnumerical-computingparallelism

Contents

Every few years someone publishes a piece titled “Fortran is dead” or “Why is Fortran still used?” — and every time, the answer is the same: because it runs things that cannot afford to be wrong. The European Centre for Medium-Range Weather Forecasts Integrated Forecasting System (IFS), which produces the best global weather predictions on Earth, is written primarily in Fortran. So is the Weather Research and Forecasting (WRF) model used operationally by national meteorological agencies worldwide. BLAS and LAPACK — the linear algebra libraries that underpin NumPy, SciPy, MATLAB, R, and nearly every ML framework that touches a dense matrix — trace their lineage directly to Fortran implementations. The DOE national laboratories run multi-decade fusion, climate, and astrophysics codes in Fortran, validated against physical experiments and peer-reviewed numerical theory. When you call numpy.dot() on a large matrix, you are, at some level, calling Fortran.

The mistake people make is conflating “old” with “obsolete.” Yes, Fortran dates to 1957. So does the transistor. The Fortran that ships in gfortran 14 or Intel ifx 2026 is not your grandfather’s FORTRAN 77. It has modules, derived types, operator overloading, object-oriented features, coarray parallelism, and a C interoperability layer precise enough to bind directly to POSIX. The FORTRAN 77 stereotype — implicit typing, column-72 source, GOTO-spaghetti, COMMON blocks — describes a dialect that has not been the standard for thirty years. The serious objection to Fortran in 2026 is not that it is ancient; it is that its ecosystem is thin, its tooling is weak, and the talent pool is shrinking. Those are real costs. But they are not the same as the language being technically inadequate for the work it does.

Why Fortran Still Owns Numerical Computing

The root cause is not inertia, though inertia plays a role. The root cause is that Fortran was designed, from the beginning, for numerical computation on arrays, and that design choice produces genuine technical advantages that still show up in modern performance measurements.

The first advantage is the absence of pointer aliasing between dummy arguments. In C, the compiler must assume that any two pointers passed to a function might point to the same memory unless you annotate them with restrict. In Fortran, dummy arguments (the Fortran term for function parameters) are defined by the standard to not alias each other — unless the programmer explicitly passes the same variable twice and the program is doing something unusual. This means the Fortran compiler can freely reorder loads and stores, vectorize loops, and keep values in registers without requiring the programmer to write restrict on every function parameter. For the inner loops of BLAS kernels and finite-difference stencils, this matters. It is one reason that Fortran BLAS implementations historically outperformed hand-ported C versions, and it remains relevant for auto-vectorization on AVX-512 and SVE targets.

The second advantage is the array language itself. Fortran arrays are first-class objects. You can write A = B + C for rank-2 arrays and the compiler generates a loop; you can write A(1:N, ::2) to slice every other column; you can pass array sections directly to subroutines. The intrinsics matmul, dot_product, sum, maxval, minval, spread, transpose, and reshape are part of the language, not a library. The compiler knows their semantics and can optimize them in ways it cannot for user-defined functions. Fortran arrays are stored in column-major order (Fortran order), which is why NumPy has an order='F' option.

The third factor is the accumulated numerical software. LAPACK, ScaLAPACK, FFTW’s Fortran interface, BLAS, LINPACK, EISPACK, NAG — these libraries represent decades of validated numerical algorithms, tested against analytical solutions and physical measurements. They are not going to be rewritten in Rust on a schedule that matters to a fusion physicist whose code has been in continuous development since 1985.

Modern Fortran: What the Standard Actually Says

The Fortran standard has evolved substantially from FORTRAN 77. The major milestones are Fortran 90/95 (free-form source, modules, dynamic allocation, array operations, interfaces), Fortran 2003 (object-oriented programming, C interoperability via iso_c_binding, IEEE arithmetic), Fortran 2008 (coarrays, do concurrent, submodules), Fortran 2018 (teams and events for coarrays, improved C interoperability with CFI_cdesc_t, further OpenMP integration), and Fortran 2023 (published November 2023), which adds conditional expressions, enumeration types, typeof() and classof(), a reduction specifier for do concurrent, improved C_F_POINTER, and removed the limit on continuation lines in free-form source.

Free-form source is the baseline for any code written after 1990. There is no column-counting; a statement can start anywhere, continue with &, and be as long as necessary. Modules replace both INCLUDE files and COMMON blocks for the sharing of data and procedure interfaces. A minimal modern module looks like this:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18


module physics_constants
  use iso_fortran_env, only: real64, int32
  implicit none
  private

  real(real64), parameter, public :: pi        = 3.141592653589793_real64
  real(real64), parameter, public :: boltzmann = 1.380649e-23_real64
  real(real64), parameter, public :: avogadro  = 6.02214076e23_real64

contains

  pure elemental function celsius_to_kelvin(t_c) result(t_k)
    real(real64), intent(in)  :: t_c
    real(real64)              :: t_k
    t_k = t_c + 273.15_real64
  end function celsius_to_kelvin

end module physics_constants

The implicit none at module scope cancels the legacy rule that variables starting with i–n are implicitly integers. That one line eliminates an entire class of latent bug. The pure attribute tells the compiler (and the reader) that this function has no side effects and does not modify any global state. The elemental attribute is particularly important: it means the function accepts either a scalar or an array of any rank, and the compiler applies it element-wise automatically. A call to celsius_to_kelvin(temperature_field) where temperature_field is a 3D array works correctly with no additional code.

Derived Types and Object-Oriented Features

Fortran 2003 introduced derived types with type-bound procedures (methods), inheritance, and polymorphism. This is not a bolt-on — the object-oriented features integrate with the array language and the pure/elemental procedure system.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32


module particle_module
  use iso_fortran_env, only: real64
  implicit none
  private

  type, public :: Particle
    real(real64) :: x, y, z
    real(real64) :: vx, vy, vz
    real(real64) :: mass
  contains
    procedure :: kinetic_energy
    procedure :: advance
  end type Particle

contains

  pure function kinetic_energy(self) result(ke)
    class(Particle), intent(in) :: self
    real(real64) :: ke
    ke = 0.5_real64 * self%mass * &
         (self%vx**2 + self%vy**2 + self%vz**2)
  end function kinetic_energy

  subroutine advance(self, dt)
    class(Particle), intent(inout) :: self
    real(real64),    intent(in)    :: dt
    self%x = self%x + self%vx * dt
    self%y = self%y + self%vy * dt
    self%z = self%z + self%vz * dt
  end subroutine advance

end module particle_module

class(Particle) rather than type(Particle) enables polymorphism: a subroutine accepting class(Particle) can receive any type that extends Particle. The percent sign (%) is the component accessor, analogous to C’s dot operator. The result is code that is clearly structured, statically typed, and still compiles to efficient machine code because the compiler can usually devirtualize calls at compile time.

The Array Language: Where Fortran Beats Naive C

Whole-array operations are more than syntactic sugar. When you write:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


! Whole-array arithmetic, no explicit loop
pressure(:,:,:) = density(:,:,:) * temperature(:,:,:) * gas_constant

! Array section: update only the interior grid points
phi(2:nx-1, 2:ny-1) = 0.25_real64 * &
    (phi(1:nx-2, 2:ny-1) + phi(3:nx,   2:ny-1) + &
     phi(2:nx-1, 1:ny-2) + phi(2:nx-1, 3:ny  ))

! Intrinsic reduction
total_energy = sum(kinetic(:) + potential(:))
max_residual = maxval(abs(residual(:,:,:)))

the compiler sees the structure explicitly and can apply SIMD vectorization, loop fusion, and cache-friendly access patterns. The same computation in C requires explicit loops that the compiler must analyze for aliasing before it can apply similar optimizations. With restrict, a good C compiler catches up; without it, or with complex pointer chains, it cannot. The Fortran version is also simply shorter and less error-prone to write.

The where construct provides array-conditional operations without a loop:

1
2
3
4


! Apply a floor to avoid negative densities (numerical diffusion artifact)
where (density < 1.0e-12_real64)
  density = 1.0e-12_real64
end where

This maps directly to a masked vector operation on hardware that supports it.

Parallelism: Coarrays, Do Concurrent, MPI, and OpenMP

Fortran’s parallelism story has multiple layers, and in practice a large HPC code typically uses several of them simultaneously alongside MPI programming essentials and OpenMP threading.

Coarray Fortran: PGAS Built Into the Language

Fortran 2008 introduced coarrays as a first-class parallel programming model using the Partitioned Global Address Space (PGAS) paradigm. A running Fortran coarray program consists of a fixed number of images — conceptually similar to MPI ranks — each with its own local memory. Coarray syntax uses square brackets to reference remote data:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24


program coarray_demo
  use iso_fortran_env, only: real64
  implicit none

  real(real64) :: local_result[*]   ! coarray: one copy per image
  real(real64) :: global_sum
  integer      :: i

  ! Each image computes its own partial result
  local_result = real(this_image(), real64) ** 2

  ! Synchronize all images before reading remote data
  sync all

  ! Image 1 gathers results from all images
  if (this_image() == 1) then
    global_sum = 0.0_real64
    do i = 1, num_images()
      global_sum = global_sum + local_result[i]
    end do
    write(*,'(a,f12.4)') 'Sum of squares: ', global_sum
  end if

end program coarray_demo

Compile and run with gfortran (OpenCoarrays backend) or ifx:

1
2
3
4
5
6
7


# gfortran with OpenCoarrays
gfortran -fcoarray=lib coarray_demo.f90 -lcaf_mpi -o coarray_demo
mpirun -np 8 ./coarray_demo

# Intel ifx (coarray support built in)
ifx -coarray coarray_demo.f90 -o coarray_demo
./coarray_demo   # or mpirun -np 8 ./coarray_demo

The PGAS model means that image 1 reading local_result[i] triggers a one-sided get from image i — no matching receive is required on image i. The communication structure is:

+----------+    +----------+    +----------+    +----------+
| Image 1  |    | Image 2  |    | Image 3  |    | Image N  |
|          |    |          |    |          |    |          |
| local_   |    | local_   |    | local_   |    | local_   |
| result   |    | result   |    | result   |    | result   |
| [1]      |    | [2]      |    | [3]      |    | [N]      |
+----+-----+    +-----+----+    +-----+----+    +-----+----+
     |                |               |               |
     +----------------+---------------+---------------+
              One-sided GET (no matching recv)
              sync all  =>  global barrier

Fortran 2018 extended coarrays with teams (subsets of images that can synchronize independently), events (non-blocking synchronization), and failed-image detection. The Fortran 2018 coarray model is more expressive than MPI in the sense that the communication is implicit in array syntax, though MPI still dominates in practice because coarray support in open-source compilers is incomplete — gfortran’s coarray support requires the OpenCoarrays library, and as of 2026 LLVM Flang does not support coarrays yet.

Do Concurrent

do concurrent was introduced in Fortran 2008 as a way to tell the compiler that loop iterations are independent and can be executed in any order:

1
2
3
4


! The compiler may vectorize, parallelize, or GPU-offload this loop
do concurrent (i = 1:n, j = 1:n)
  c(i,j) = a(i,j) * b(i,j) + alpha
end do

Fortran 2023 adds a reduce clause to do concurrent, closing a significant gap with OpenMP’s reduction directive. Intel ifx and NVIDIA nvfortran can both target GPU execution from do concurrent — ifx maps it to Intel GPU kernels and nvfortran maps it to CUDA through the NVIDIA HPC SDK. This is meaningful: it means a loop annotated once can be offloaded without rewriting it as an OpenACC or OpenMP target region. See GPU programming without CUDA for the broader GPU offload landscape.

Relationship to MPI and OpenMP

Real HPC codes use a hybrid: MPI for distributed-memory parallelism across nodes, OpenMP for shared-memory threading within a node, and increasingly GPU offload directives for accelerators. Fortran is a full citizen in all three models. The MPI Fortran interface is part of the MPI standard; OpenMP Fortran directives (!$omp parallel do, !$omp simd) work identically to their C counterparts; OpenACC is primarily a Fortran/C feature. Running hybrid MPI+OpenMP Fortran jobs on a cluster means the SLURM job scheduler setup is identical to what you’d use for C or C++ — --ntasks-per-node, --cpus-per-task, and OMP_NUM_THREADS are language-agnostic.

C Interoperability: iso_c_binding

Since Fortran 2003, the standard intrinsic module iso_c_binding provides a precise, portable ABI bridge between Fortran and C. Before this existed, calling Fortran from C or vice versa was compiler-specific and fragile. Now it is well-defined.

Calling a C function from Fortran:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


module c_interfaces
  use iso_c_binding
  implicit none

  interface
    ! Binds to C: int posix_memalign(void **memptr, size_t align, size_t size)
    function posix_memalign(memptr, alignment, sz) &
        bind(c, name='posix_memalign') result(rc)
      import :: c_ptr, c_size_t, c_int
      type(c_ptr),    intent(out) :: memptr
      integer(c_size_t), value   :: alignment
      integer(c_size_t), value   :: sz
      integer(c_int)             :: rc
    end function posix_memalign
  end interface

end module c_interfaces

Exposing a Fortran subroutine to C:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12


! This subroutine has C linkage and can be called from C as:
! void dgemm_wrapper(double *a, double *b, double *c, int n);
subroutine dgemm_wrapper(a, b, c, n) bind(c, name='dgemm_wrapper')
  use iso_c_binding
  implicit none
  integer(c_int), value        :: n
  real(c_double), intent(in)   :: a(n,n), b(n,n)
  real(c_double), intent(inout):: c(n,n)

  ! Call BLAS DGEMM
  call dgemm('N','N', n, n, n, 1.0d0, a, n, b, n, 0.0d0, c, n)
end subroutine dgemm_wrapper

Key points: bind(c) enforces C calling convention (value semantics for scalars unless intent implies otherwise, C name mangling). The value attribute on scalar arguments passes by value rather than by reference, matching C convention. Kind parameters from iso_c_binding — c_int, c_double, c_size_t, c_ptr — map to the corresponding C types on any conforming platform.

Fortran 2018 extended this with CFI_cdesc_t, a C descriptor for Fortran arrays, allowing C code to manipulate Fortran assumed-shape arrays without knowing their exact memory layout at compile time. This is used by f2py and similar tools under the hood.

Calling Fortran from Python

The three main approaches are f2py, ctypes, and cffi.

f2py (part of NumPy) remains the standard tool for wrapping Fortran subroutines and functions for Python consumption. It works well for F77-style and F90 code, generates a C extension module, and handles type mapping automatically for numeric arrays. Its limitations are real: derived types are not fully supported, assumed-shape arrays require extra annotation, and modules with complex use-association can confuse the parser.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


# Simple f2py workflow for a Fortran module
# Given a file 'numerical.f90' with subroutines:

# Step 1: generate a signature file to inspect/annotate
python -m numpy.f2py numerical.f90 -m numerical_ext --overwrite-signature

# Step 2: compile to a Python extension
python -m numpy.f2py -c numerical.f90 -m numerical_ext \
    --f90flags="-O3 -march=native"

# Step 3: use from Python
# import numerical_ext
# result = numerical_ext.my_subroutine(array_arg)

With gfortran and Meson as the f2py build backend (available since NumPy 1.22):

1
2


python -m numpy.f2py -c numerical.f90 -m numerical_ext \
    --backend meson --f90flags="-O3"

ctypes is the fallback when f2py fails. Compile the Fortran code with bind(c) interfaces into a shared library, then load it with ctypes.CDLL. The iso_c_binding machinery makes this precise:

1

gfortran -shared -fPIC -O3 my_routines.f90 -o libmy_routines.so

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


import ctypes, numpy as np

lib = ctypes.CDLL('./libmy_routines.so')
lib.my_kernel.argtypes = [
    ctypes.POINTER(ctypes.c_double),
    ctypes.POINTER(ctypes.c_double),
    ctypes.c_int
]
lib.my_kernel.restype = None

a = np.asfortranarray(np.ones((100,100)))
b = np.asfortranarray(np.zeros((100,100)))
lib.my_kernel(a.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
              b.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
              ctypes.c_int(100))

Note np.asfortranarray(): Fortran column-major ordering must match. Passing a C-order array to a Fortran routine that assumes column-major layout is a silent correctness bug that produces wrong numbers without any error.

The Julia scientific computing ecosystem provides a complementary approach: Julia’s ccall can invoke Fortran functions directly (with bind(c) interfaces), and Julia’s type system maps naturally to Fortran’s numeric kinds. For new code that needs to call both legacy Fortran libraries and leverage Julia’s JIT performance, this is increasingly attractive.

The Compiler Landscape in 2026

Compiler	License	F2018	F2023	GPU Offload	Coarrays	Notes
gfortran (GCC 14+)	GPL	Full	Partial	None (no built-in)	Via OpenCoarrays	The universal workhorse; best for portability and CI
LLVM Flang (LLVM 20/21)	Apache 2.0	Mostly complete	Early	Via OpenMP target	Not yet	Renamed from flang-new in LLVM 20 (March 2025); OpenMP now stable; coarrays and PDTs still missing
Intel ifx (oneAPI 2026)	Proprietary	Full	Several features	Intel GPU (OpenMP target, do concurrent)	Full	ifort EOL’d in 2025; ifx is the replacement; excellent on Xeon; MKL integration
NVIDIA nvfortran (HPC SDK 26.x)	Proprietary	Full	Partial	NVIDIA GPU (OpenACC, OpenMP target)	Limited	Best-in-class for NVIDIA GPU offload; OpenACC 2.7
LFortran	BSD	Partial	Early	Experimental	No	Interactive/REPL use; near-beta in early 2026 after compiling fpm
Cray CFE (HPE)	Proprietary	Full	Partial	Via CCE/OpenMP	Full	Used on Cray/HPE Slingshot clusters; strong coarray implementation

The LLVM Flang story deserves elaboration. For years, the official LLVM Fortran compiler was called flang-new to distinguish it from an older, unrelated Flang project. LLVM 20, released in March 2025, dropped the -new suffix: the binary is now simply flang. OpenMP support graduated from experimental to stable in that release cycle. The two significant missing pieces as of mid-2026 are coarrays and parameterized derived types with length parameters. These are not trivial omissions — coarrays require runtime support, and PDTs have complex type-system interactions — but for the large class of codes that use neither, Flang is genuinely usable and benefits from the entire LLVM optimization pipeline, including polyhedral loop optimizations and the same vectorization backends that power Clang.

Compiling the same code with different compilers:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14


# gfortran — most portable
gfortran -O3 -march=native -ffast-math physics.f90 -o physics_gfortran

# LLVM Flang
flang -O3 -march=native physics.f90 -o physics_flang

# Intel ifx with MKL
ifx -O3 -xHost -qmkl physics.f90 -o physics_ifx

# NVIDIA nvfortran with GPU offload via OpenACC
nvfortran -O3 -acc=gpu -gpu=cc90 physics.f90 -o physics_nvfortran

# Intel ifx with do concurrent GPU offload
ifx -O3 -xHost -fopenmp-target-do-concurrent physics.f90 -o physics_ifx_gpu

Reading and Modernizing Legacy Fixed-Form Code

A large fraction of production Fortran code in use today is still in fixed-form FORTRAN 77 or mixed fixed/free F90. Understanding it is a practical necessity.

Fixed-form has strict column rules: columns 1–5 are for statement labels, column 6 is the continuation marker (any non-blank, non-zero character), columns 7–72 are the statement, and columns 73–80 are traditionally ignored (punch-card sequence numbers). A comment line starts with C or * in column 1.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


C     This is a fixed-form comment (column 1)
      PROGRAM LEGACY
      IMPLICIT INTEGER (I-N)       ! default implicit typing still active
      REAL A(100), B(100)
      COMMON /SHARED/ A, B         ! global data via COMMON block
      DO 100 I=1,100
        A(I) = REAL(I) * 3.14159
  100 CONTINUE
      GOTO 200
      WRITE(*,*) 'unreachable'
  200 CALL MYSUB(A, B, 100)
      STOP
      END

When modernizing incrementally, the highest-leverage changes are:

Add IMPLICIT NONE to every program unit immediately. This causes the compiler to flag all previously-implicit variable declarations; the resulting errors are a map of bugs.
Convert COMMON blocks to modules. Data in a COMMON block has no type checking at the boundaries — any program unit can declare it differently and get away with it until runtime. A module variable is checked everywhere it is used.
Replace labeled DO/CONTINUE loops with modern do/end do.
Replace GOTO with structured control flow. Most FORTRAN 77 GOTO usage implements what are now if/else if, exit, cycle, or select case.
Add explicit interfaces (or move procedures into modules). Without an explicit interface, the compiler cannot check argument types or array shapes at call sites.

You do not have to do this all at once. A COMMON block can be replaced by a module, and the rest of the code is left unchanged, because a Fortran program can mix fixed-form and free-form source files in the same link step (though not in the same file).

The Honest Assessment: Fortran’s Weaknesses

The ecosystem is genuinely thin. The Fortran Package Manager (fpm) exists and released version 0.12.0 in May 2025 with compile_commands.json export, BLAS/LAPACK metapackages, and HDF5 support. It is well-designed and actively developed. But the registry of packages is nothing like what Cargo or pip offers. If you need string parsing, a JSON library, an HTTP client, or a GUI, you are writing FFI bindings or maintaining a C dependency. These are not hard problems, but they are additional work that does not exist in Python or even in Julia’s ecosystem.

Fortran’s string handling is painful. Character variables have fixed declared length, or use allocatable character arrays, or use the varying_string type that almost nobody uses. There is no standard regular expression support. File I/O with formatted records works fine for scientific data, but anything involving variable-format text, CSV parsing, or structured data (JSON, YAML, HDF5 via direct Fortran) requires either a third-party library or dropping out to C.

Error handling is weak. Fortran has no exceptions. The standard STAT= and ERRMSG= specifiers on allocation and I/O provide a mechanism for checking errors, but propagating them through a call stack requires passing status arguments by hand — a pattern that disappears from codebases because it is tedious, leaving unhandled error paths throughout real scientific codes.

The tooling gap is significant. There is no Fortran language server with the maturity of clangd or rust-analyzer, though fortls (the Fortran Language Server) is usable and improving. Debugger integration with modern Fortran features — particularly allocatable arrays and derived types — is uneven in gdb. Static analysis tools are sparse. IDE support is second-class everywhere except with Intel’s Visual Studio integration.

The talent pool is contracting. Universities have largely stopped teaching Fortran. The people who can read a 40,000-line fixed-form climate model and understand it are disproportionately approaching retirement. New HPC codes are increasingly started in C++, and sometimes in Julia. This does not make Fortran wrong for the job — it makes organizational risk management a legitimate concern when choosing it for new projects.

Verdict

Fortran is not dead. It is not dying quickly. The world’s most consequential numerical simulations run on it, and they will not be rewritten this decade or the next. The language standard is genuinely modern: iso_c_binding, modules, derived types, do concurrent, coarrays, and elemental procedures are not legacy features — they are language-level tools for writing correct, fast, parallel numerical code. The compiler ecosystem in 2026 is the strongest it has been in twenty years: gfortran is mature and ubiquitous, LLVM Flang is production-viable for the large class of codes that do not use coarrays or PDTs, Intel ifx replaces ifort cleanly with better standards conformance and GPU offload, and NVIDIA nvfortran handles the CUDA-offload use case better than any alternative.

The honest constraints are real: you are signing up for a narrow ecosystem, weak string ergonomics, non-existent package infrastructure by modern standards, and a shrinking pool of developers who can maintain the code. For new numerical projects that do not carry legacy dependencies, the cost-benefit calculation against C++ or Julia deserves honest evaluation.

For the codes that already exist — the IFS, WRF, LAPACK, the national laboratory multiphysics codes — the question is not “should we use Fortran?” but “how do we keep this running correctly on the next generation of hardware?” The answer involves ifx for Intel GPU offload, nvfortran for NVIDIA, do concurrent where coarrays are overkill, and iso_c_binding to expose validated numerical kernels to the Python data science stack above. That is not archaeology. That is engineering.