Blog › ICP guides

Fortran developer on retainer: HPC, OpenMP, MPI, BLAS/LAPACK, and scientific computing on monthly retainer

September 17, 2026 · ~20 min read

A climate modeling group had been running their global atmospheric simulation on gfortran for four years. When they migrated to ifort on a new cluster to take advantage of its AVX-512 code generation, the model passed all physical validation checks but failed to reproduce bit-identical results against the reference gfortran output. This mattered: the publication record for the model cited specific numerical checksums from ensemble runs, and reviewers expected any code change that altered floating-point results to come with an explicit reproducibility statement. The team had changed no physics, no algorithm, and no numerical method. The only change was the compiler.

The Fortran developer on retainer diagnosed it in two separate passes over three days. The first pass compared generated assembly with -S: gfortran with default flags emitted VADDSD at every multiply-accumulate site in the inner loops; ifort emitted VFMADD231PD — fused multiply-add, with the intermediate product unrounded. Adding -ffp-contract=off to the ifort flags suppressed FMA contraction and brought results closer, but not bit-identical. The second pass examined the COMMON block structures shared between compilation units. A GlobalState COMMON block in atm_physics.f90 lacked the SEQUENCE attribute while climate_dynamics.f90 used it without declaring it — without SEQUENCE, field alignment and padding in a derived type inside a COMMON block are compiler-defined; the two compilers had placed a four-byte padding field in different positions, causing different bytes to be read as the next variable in a loop over the block's memory. Adding SEQUENCE to GlobalState in both compilation units, combined with -ffp-contract=off, produced bit-identical results. Total work: approximately 22 hours across three sessions of assembly comparison and COMMON block auditing.

No physics changed. No algorithm changed. Two compiler flags and one SEQUENCE attribute keyword added across three source files. The invisible artifact was the elimination of a floating-point reproducibility failure that would have prevented peer-reviewed publication of the model's next ensemble run. A Fortran developer on monthly retainer does this category of work continuously: auditing floating-point contraction settings before compiler migrations break reproducibility, optimizing LAPACK workspace allocation before diagonalization throughput limits ensemble throughput, and diagnosing OpenMP thread oversubscription before it degrades HPC node efficiency to below-serial performance.

Modern Fortran: allocatables, modules, assumed-shape arrays, and KIND portability

The transition from legacy Fortran 77 idioms to modern Fortran 90 through 2018 features is the most common architectural work in a Fortran retainer. COMMON blocks are the primary target: they provide no type checking, no scoping, and no protection against field-order inconsistency across compilation units. Modules with USE association replace them with explicit, type-checked imports. The USE physics_constants, ONLY: gravity, earth_radius pattern imports only what is needed from the module, making dependencies visible in the source rather than implicit in the COMMON block layout. Module variables can carry the PROTECTED attribute, making them readable from any USE-associating scope but writeable only from within the defining module — enforcing a read-only interface for constants and reference values that should not be modified by downstream code. Interface blocks defined in a module give the compiler full argument type and shape information for every procedure call, enabling argument mismatch errors at compile time rather than silent wrong-value bugs at runtime.

Allocatable arrays are the modern replacement for fixed-size arrays dimensioned at compile time or passed as explicit-shape dummy arguments. The lifecycle — ALLOCATE(a(n,m), stat=ierr) with ierr checked for non-zero, operations on a, DEALLOCATE(a) — is straightforward, but retainer work concentrates on the cases around it: checking ALLOCATED(a) before an allocation that may be called multiple times prevents double-allocation errors; automatic deallocation at scope exit for local allocatable variables (introduced in Fortran 95) means that allocatables declared inside a subroutine are freed when the subroutine returns, eliminating the memory leak class common in manually managed explicit-shape arrays; and reallocatable assignment in Fortran 2003 — where a = b for two allocatable arrays of different shape automatically reallocates a to match the shape of b before copying — enables clean array resizing without explicit DEALLOCATE/ALLOCATE sequences. Assumed-shape dummy arguments, declared as real(8), intent(in) :: a(:,:), eliminate the need to pass explicit dimensions as separate arguments; the shape is carried with the array descriptor. The CONTIGUOUS attribute on an assumed-shape or pointer dummy argument tells the compiler that the array occupies contiguous memory — enabling alias analysis and making the array eligible for SIMD vectorization that the compiler cannot apply to a non-contiguous assumed-shape argument.

KIND parameter portability is the least visible and most consequential design decision in a Fortran codebase that runs on multiple platforms. Hard-coding real*8 or real(8) as the double precision kind is not portable across compilers and architectures — the correct approach is integer, parameter :: dp = selected_real_kind(15,307), which requests the smallest KIND that provides at least 15 decimal digits of precision and a decimal exponent range of at least 307, and then using real(dp) throughout. Propagating this KIND parameter consistently across a large codebase requires a retainer audit: every literal constant that feeds a double-precision computation must carry the kind suffix (1.0_dp rather than 1.0 or 1.0d0), every external procedure that returns a double-precision value must have an explicit interface that declares the return kind, and every intrinsic function call must use the KIND-tagged argument to ensure the intrinsic operates at the right precision. Elemental and pure function design reinforces this: an elemental function declared elemental real(dp) function scale(x, factor) operates element-wise across any rank array without explicit loops, and its pure attribute (no side effects, no I/O, no modification of non-local state) enables the compiler to parallelize its application across array elements and evaluate it in the context of do concurrent loops.

BLAS/LAPACK, FFTW3, NetCDF/HDF5, and numerical library integration

BLAS and LAPACK are the numerical backbone of most Fortran HPC codebases, and the retainer work around them concentrates on three areas: correct external declarations, argument order discipline, and workspace optimization. External declarations for BLAS/LAPACK routines — external dgemm, dgesv, dsyev or a module interface that provides the full signature — are necessary because these routines are Fortran 77-era procedures without modern interface blocks. The dgemm call for C = alpha*A*B + beta*C takes 13 arguments in order: TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC — transposing TRANSA and TRANSB silently computes a different product without an error. The dgesv call for AX = B returns an INFO code that is 0 for success, negative for an illegal argument (with |INFO| giving the argument position), and positive for a singular U factor (with INFO giving the diagonal index where it occurred) — codebases that do not check INFO after dgesv will silently use the solution of a singular system. For dsyev, the correct pattern requires two calls: the first with LWORK=-1, which writes the optimal workspace size as a real number into WORK(1) without computing the eigenvalues, allowing the caller to ALLOCATE(WORK(int(WORK(1)))) before the second call that actually computes the eigenvalues. Codebases that skip the query call and use a fixed LWORK — often chosen conservatively as 3*N — leave performance on the table: the optimal workspace for large symmetric matrices is typically 6*N to 8*N, and the underprovided workspace degrades the blocked algorithm to a less efficient unblocked path.

FFTW3's Fortran interface uses the ISO C binding: use iso_c_binding and include 'fftw3.f03' provide the Fortran-compatible type definitions and subroutine interfaces. A one-dimensional real-to-complex transform requires fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE) to create a plan, fftw_execute_dft_r2c(plan, in, out) to execute it, and fftw_destroy_plan(plan) to release the plan memory. The planning flag FFTW_ESTIMATE makes a quick heuristic plan; FFTW_MEASURE benchmarks candidate algorithms and selects the fastest, which can be 2–5x faster for large transforms but takes seconds at startup. Wisdom import and export — fftw_import_wisdom_from_filename/fftw_export_wisdom_to_filename — serializes the measurement results to a file so that subsequent runs skip the measurement phase while using the measured-optimal plan. For OpenMP-parallel FFT, fftw_plan_with_nthreads(omp_get_max_threads()) before the plan creation enables multi-threaded execution; the plan is thread-safe for execution (multiple threads may call fftw_execute_dft_r2c concurrently on the same plan with different in/out buffers), but plan creation is not thread-safe and must occur in a serial region.

NetCDF and HDF5 are the standard storage formats for climate, atmospheric, and oceanographic model output. The NetCDF Fortran 90 API — use netcdf — provides nf90_open/nf90_create/nf90_close for file lifecycle, nf90_inq_varid and nf90_inq_dimid for variable and dimension lookup by name, and nf90_get_var/nf90_put_var with start and count integer arrays for hyperslab access — reading a time slice from a four-dimensional (lon, lat, level, time) variable uses start=[1,1,1,t] and count=[nlon,nlat,nlev,1]. An unlimited time dimension enables append-writes across ensemble members without pre-sizing the time axis. Defining variables with nf90_def_var(..., deflate_level=4) enables zlib compression at level 4, which typically reduces output size by 30–60% for floating-point climate fields with moderate spatial correlation. HDF5's Fortran bindings require explicit initialization and finalization: h5open_f(ierr) and h5close_f(ierr) bracket the entire HDF5 usage in a program; each file, dataset, and dataspace has its own integer(hid_t) handle that must be explicitly closed with the matching h5fclose_f, h5dclose_f, or h5sclose_f to avoid handle leaks that accumulate across long simulation runs. A retainer engagement covering HDF5 integration audits every hid_t handle for a matching close call and verifies that the file is closed cleanly under both normal-completion and error-exit code paths.

OpenMP parallelism, MPI distribution, and compiler-driven performance optimization

OpenMP in Fortran HPC codebases is deceptively easy to introduce and expensive to get right. The !$omp parallel do directive with schedule(static) distributes loop iterations in equal-sized blocks across threads — correct for uniform-cost iterations like simple array arithmetic but wrong for loops where iteration cost varies, where schedule(dynamic,chunk) assigns a new block of chunk iterations to each thread as it finishes its previous block. The data-sharing clauses require explicit thought: firstprivate(x) gives each thread a private copy of x initialized to the value of x in the enclosing scope; lastprivate(i) writes the value of i from the last loop iteration back to the enclosing scope when the parallel region exits; reduction(+:sum) gives each thread a private zero-initialized sum, accumulates locally, and combines all thread-local sums into the shared sum after the loop. The nowait clause on a worksharing construct removes the implicit barrier at the end of that construct, allowing threads to proceed to the next worksharing construct without synchronizing — useful when two independent array operations follow each other within a parallel region. The OpenBLAS nested-thread oversubscription hazard is the most common performance regression in OpenMP-enabled HPC Fortran: if OpenBLAS was compiled with OpenMP support and the Fortran program is already running inside a !$omp parallel region, each BLAS call spawns OPENBLAS_NUM_THREADS additional threads, multiplying the thread count by the BLAS thread count. Setting OPENBLAS_NUM_THREADS=1 (or calling openblas_set_num_threads(1) from Fortran via the C binding) before entering the parallel region limits each BLAS call to a single thread within the OpenMP worker threads, eliminating the oversubscription.

MPI domain decomposition is the primary mechanism for distributing Fortran simulations across multiple nodes. The baseline setup — call mpi_init(ierr), call mpi_comm_rank(MPI_COMM_WORLD, rank, ierr), call mpi_comm_size(MPI_COMM_WORLD, nprocs, ierr), call mpi_finalize(ierr) — is the scaffold around which domain decomposition logic is built. Blocking mpi_send/mpi_recv is the simplest communication pattern but serializes computation and communication: the sending rank blocks at mpi_send until the receiving rank posts its mpi_recv. Non-blocking mpi_isend/mpi_irecv return immediately with a MPI_Request handle; the program continues computing while the message is in transit, and call mpi_waitall(count, requests, MPI_STATUSES_IGNORE, ierr) blocks only when the message data is actually needed. This overlap of communication and computation is the primary performance lever in latency-bound domain decomposition. Global reductions — summing a scalar across all ranks — use call mpi_allreduce(local_sum, global_sum, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ierr); mpi_allreduce makes the result available on every rank, while mpi_reduce makes it available only on rank 0. Sending a Fortran derived type with mixed field types (integer and real fields in the same structure) requires a derived datatype: mpi_type_create_struct takes arrays of block lengths, byte displacements (computed with mpi_get_address), and old MPI datatypes, producing a new MPI datatype that describes the layout of the derived type in memory; the type must be committed with mpi_type_commit before use and freed with mpi_type_free when no longer needed.

Compiler-driven performance optimization in Fortran follows a three-stage process that retainer work enforces systematically. First, establish the baseline: -O3 -march=native enables aggressive optimization and generates instructions tuned to the native CPU. Second, diagnose vectorization: gfortran -fopt-info-vec-optimized -fopt-info-vec-missed reports which loops were vectorized and which were not and why; ifort's -qopt-report=5 produces a detailed per-loop report including the reason for non-vectorization (aliasing assumed, non-contiguous memory access, loop-carried dependency detected). The most common vectorization blockers in Fortran are assumed-shape pointer arguments without the CONTIGUOUS attribute (the compiler assumes possible non-contiguous strides), stride-1 access violations from column-major array traversal in row-major order (Fortran arrays are column-major; the innermost loop index must be the first array index for stride-1 access and cache efficiency), and loop-carried dependencies introduced by shared reduction variables that should be expressed with reduction clauses. Third, apply profile-guided optimization: compile with -fprofile-generate, run a representative workload to produce .gcda profile data files, then recompile with -fprofile-use to apply the measured branch probabilities and hot-loop information to the optimization decisions. For floating-point contract consistency across compilers, -ffp-contract=fast enables FMA fusion (approximately 2–5% throughput gain for dense linear algebra) at the cost of cross-compiler non-reproducibility; -ffp-contract=off suppresses FMA contraction and produces bit-identical results between gfortran and ifort for the same source code, which is the correct setting for any codebase where reproducibility is a publication or validation requirement.

How HourTab tracks Fortran developer retainer hours

Fortran developer retainers produce some of the most disproportionate work-to-visible-output ratios in any language retainer. A session that diagnosed a floating-point reproducibility failure between gfortran and ifort produced two flag changes and one SEQUENCE keyword across three source files. The session involved compiling both toolchains with -S to compare generated assembly, identifying the VFMADD231PD vs VADDSD FMA contraction discrepancy, adding -ffp-contract=off to the ifort flags and observing residual non-identity, auditing COMMON block structures for SEQUENCE attribute consistency across compilation units, and confirming bit-identical results against the reference checksums. The log entry “fixed compiler reproducibility, 22h” gives the client no path from 22 hours to the elimination of the reproducibility failure that blocked ensemble publication — because nothing in two flag changes and one keyword communicates the assembly comparison and COMMON block audit that produced them.

HourTab gives Fortran developers a public retainer-hours URL they share with each client at the start of the engagement. The client opens the URL and sees the current burn-down without logging in. For Fortran retainers specifically, the work log format carries the weight: each entry should name the LAPACK routine and its INFO return code when relevant (dsyev INFO=0 after fix; INFO=3 — singular U factor at index 3 — before the matrix conditioning improvement), the compiler flag that was incorrect and what it generated (ifort default: -ffp-contract=fast, VFMADD231PD; after fix: -ffp-contract=off, VADDSD — bit-identical with gfortran), the OpenMP clause that was missing or wrong (reduction(+:sum) added to inner loop — was accumulating into shared variable, data race; -fopt-info-vec confirms vectorization enabled after fix), the MPI pattern changed and the overlap gain (blocking mpi_send/mpi_recv replaced with mpi_isend/mpi_irecv + mpi_waitall; computation-communication overlap: 34% wall time reduction on 128-rank job), the scope of the audit, and the before/after metric. That entry takes five minutes to write and turns the client check-in from a twenty-minute explanation of what FMA contraction is and why SEQUENCE matters for COMMON block portability into a two-sentence acknowledgment of the reproducibility confirmation.

The retainer model fits Fortran HPC platform engineering because the numerical, compiler, and parallelism landscape continues to evolve — new HPC clusters bring new CPU microarchitectures requiring -march=native re-validation, new gfortran and ifort major releases change default optimization behavior, and new MPI implementations alter collective algorithm selection that affects domain decomposition performance. A project contract closes when the current reproducibility audit or LAPACK optimization milestone is complete. A Fortran retainer stays open for the next compiler upgrade that changes FMA behavior, the next OpenMP addition that introduces a scheduling imbalance, and the next MPI domain decomposition redesign triggered by a 2x increase in ensemble size.

FAQ: Fortran developer retainers