Blog › ICP guides

Fortran developer on retainer: array indexing bounds, assumed-shape arrays, OpenMP parallelization, and Fortran on monthly retainer

October 2, 2026 · ~14 min read

A Fortran developer was maintaining an atmospheric climate model. The model included a radiation subroutine that computed temperature profiles through vertical layers of the atmosphere. After a refactoring pass that extracted the vertical indexing into a separate array, one simulation output showed a vertical temperature profile that was shifted one level: the surface temperature appeared at level two, and the top-of-atmosphere temperature was missing. The other 11 grid-point outputs in that batch were correct.

The main program declared the temperature array with an explicit 0-based lower bound:

REAL :: temperature(0:nz-1, ny, nx)

The radiation subroutine accepted the array as an assumed-shape dummy argument:

SUBROUTINE compute_radiation(temperature, ...)
  REAL, INTENT(IN) :: temperature(:,:,:)

In Fortran, an assumed-shape dummy argument ((:,:,:)) receives the actual shape and bounds from the calling unit. When the main program’s temperature(0:nz-1, ny, nx) array was passed, the subroutine’s temperature dummy argument had bounds (0:nz-1, 1:ny, 1:nx) — inheriting the explicit zero lower bound in the first dimension. The DO loop inside the subroutine iterated from 1 to SIZE(temperature, 1). SIZE(temperature, 1) returned nz (the number of elements), but the loop starting at index 1 began at the second element of the 0-indexed first dimension, skipping level 0 (the surface) and reading through level nz-1 (one past the intended top). The surface temperature was never read; the top level was read from out-of-declared-bounds memory. One simulation output with a shifted vertical profile. The fix: change the main-program declaration to use default 1-based indexing — REAL :: temperature(nz, ny, nx) — and update the allocation and access sites to use 1-based indices consistently. Wrong simulation outputs: 1 → 0. The work log said “fixed vertical profile shift in radiation, 6h.” What the log does not say is that Fortran assumed-shape arrays inherit the calling unit’s lower bound — and that mixing explicit 0-based lower bounds in the calling unit with implicit 1-based loops in subroutines is the most common source of off-by-one errors in Fortran refactoring.

Fortran overview: the language of numerical simulation since 1957

Fortran (Formula Translation) was developed at IBM by John Backus and his team and released in 1957 as the first high-level programming language to achieve widespread adoption. The original goal was to produce code that ran as fast as hand-written assembly while being expressible in mathematical notation. The first Fortran compiler generated code that matched or exceeded assembly programmer output for numerical computations, establishing Fortran as the language of scientific computing for decades. Fortran 77 standardized fixed-format source code and added structured control flow. Fortran 90 added free-format source, modules, derived types, and the modern array syntax. Fortran 95 refined array operations. Fortran 2003 added object-oriented features, interoperability with C, and improved I/O. Fortran 2008 added DO CONCURRENT for data-parallel loops and coarrays for distributed-memory parallel programming. Fortran 2018 is the current standard, adding further parallel and C interoperability features.

Fortran remains the dominant language for computational science because the codebase is irreplaceable. Climate models, weather forecasting models, plasma physics simulations, fluid dynamics solvers, structural engineering finite element codes, and quantum chemistry programs have been developed and validated over decades in Fortran. The validation history — the knowledge that a particular Fortran subroutine produces results that match physical measurements — is as important as the code itself. Rewriting a validated Fortran climate model in Python or Julia requires re-validating every subroutine from scratch, which can take longer than the original development. The models run. The science depends on them running correctly. The maintenance work falls on Fortran developers who understand the language’s array model and numerical behavior.

Fortran retainer work today covers: climate and weather model maintenance (NCAR’s WRF, GFDL models, ECMWF IFS — maintained by national laboratories and weather forecasting agencies); computational chemistry maintenance (Gaussian, GAMESS, NWChem — maintained by university chemistry departments and pharmaceutical companies); finite element analysis code maintenance (NASTRAN, ABAQUS user subroutines, OpenFOAM Fortran components); HPC performance optimization (tuning Fortran codes for specific cluster architectures, GPU acceleration via CUDA Fortran or OpenACC, OpenMP and MPI parallelization); and compiler migration (moving from legacy fixed-format Fortran 77 or Fortran II code to modern free-format Fortran 90/2018 with modules and explicit interfaces).

Fortran array model: 1-based indexing, explicit bounds, and assumed-shape arguments

Fortran arrays are 1-indexed by default. An array declared as REAL :: pressure(n) has valid indices 1 through n. An array declared as REAL :: pressure(10, 20) has valid indices 1 through 10 in the first dimension and 1 through 20 in the second. Fortran stores multi-dimensional arrays in column-major order (the first index varies fastest in memory), which is the opposite of C’s row-major order. A Fortran developer working with a two-dimensional array accesses it as A(row, col), but in memory the elements are A(1,1), A(2,1), A(3,1), ... — columns stored consecutively.

Explicit lower bounds are allowed: REAL :: altitude(0:nz-1) declares an array with valid indices 0 through nz-1. Explicit lower bounds are common in code written by developers with C backgrounds, or in code that represents data with a natural 0-based indexing (e.g., a frequency bin array where bin 0 is the DC component). The bug: explicit lower bounds in the calling unit interact with subroutine dummy argument declarations in non-obvious ways depending on whether the dummy argument is explicit-shape, assumed-shape, or assumed-size.

An explicit-shape dummy argument (REAL :: arr(n)) has bounds 1:n regardless of the calling unit’s bounds. A calling unit that passes temperature(0:nz-1) to an explicit-shape dummy temperature(nz) will have the subroutine see indices 1:nz, where index 1 corresponds to the calling unit’s index 0. This is a well-known idiom in Fortran 77 code, but produces index-shifting bugs when developers assume the subroutine’s index 1 corresponds to the calling unit’s index 1.

An assumed-shape dummy argument (REAL :: arr(:)) inherits the bounds from the calling unit. If the calling unit passes temperature(0:nz-1), the assumed-shape dummy argument has bounds 0:nz-1 in the subroutine. This is the modern Fortran 90 convention and is generally correct — but only if the subroutine’s loops use the dummy argument’s LBOUND and UBOUND intrinsics to determine loop limits, rather than hardcoding start at 1. A loop that starts at 1 when the assumed-shape argument has lower bound 0 skips the first element. The diagnosis: check the lower bound of the actual argument in the calling unit; check whether the subroutine’s loops use LBOUND(arr,1) or hardcode 1 as the start.

An assumed-size dummy argument (REAL :: arr(*)) is the Fortran 77 style; the subroutine cannot determine the array size and must be told via a separate scalar argument. Assumed-size arguments cannot be used with Fortran 90 array intrinsics like SIZE. Code that mixes assumed-size dummy arguments with calls to SIZE will not compile or will produce wrong size values at runtime depending on the compiler’s handling.

OpenMP parallelization, race conditions, and IMPLICIT NONE discipline

OpenMP adds shared-memory parallelism to Fortran via compiler directives. A parallel loop is annotated with !$OMP PARALLEL DO; multiple threads execute the loop body concurrently. The most common Fortran OpenMP bug: a loop accumulator variable that is semantically private to each iteration is declared as SHARED (the default for variables that appear in the loop but are not loop variables). Multiple threads read and write the shared accumulator simultaneously without synchronization, producing a non-deterministic result that varies between runs.

The correct fix depends on the accumulator’s role. If each iteration computes an independent accumulator value (e.g., a running total for that iteration only), the accumulator should be declared PRIVATE: !$OMP PARALLEL DO PRIVATE(local_sum). If the accumulator is a grand total that all iterations contribute to (e.g., computing a global sum over all array elements), it should be declared with a REDUCTION clause: !$OMP PARALLEL DO REDUCTION(+:total). The REDUCTION clause creates a private copy of the variable in each thread, accumulates each thread’s partial result, and combines them correctly after the parallel region. A developer who adds OpenMP to an existing serial loop without analyzing which variables are per-iteration (PRIVATE) vs cross-iteration (REDUCTION) will produce wrong results whenever multiple threads update the same accumulator concurrently.

IMPLICIT NONE is the Fortran declaration that disables implicit typing — the historical Fortran rule that variables beginning with letters I through N are INTEGER and all others are REAL, without any declaration. Legacy Fortran code without IMPLICIT NONE can have typos in variable names that silently create new implicitly-typed variables. A developer who writes temperture instead of temperature in a loop body does not get a compile error in code without IMPLICIT NONE; the compiler creates a new variable temperture initialized to zero. The loop writes to the new variable; the intended array is never updated. Adding IMPLICIT NONE to all program units is a retainer work task that converts these silent bugs into compile-time errors — typically surfacing 3 to 8 existing variable name typos per legacy subroutine.

BLAS and LAPACK interfaces, column-major layout, and leading dimension errors

BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra PACKage) are the standard numerical linear algebra libraries for Fortran and scientific computing. BLAS Level 3 routines like DGEMM (double-precision general matrix multiply) are the core of high-performance linear algebra. The BLAS Fortran API uses column-major matrix layout (matching Fortran’s native array order) with a leading dimension parameter (LDA) that specifies the number of rows in the full array as allocated, allowing submatrix views to be passed without copying.

The most common BLAS retainer bug: a developer passes a submatrix to DGEMM but provides the submatrix row count as LDA instead of the full array’s row count. A matrix declared as REAL*8 :: A(1000, 800) with only the top-left 300 × 250 submatrix used for a particular call must be passed with LDA = 1000 (the declared row dimension), not LDA = 300 (the submatrix row count). With the wrong LDA, BLAS reads column offsets incorrectly — each column starts 700 rows too early in memory — producing a matrix multiply result that is numerically wrong with no error message. The diagnostic: check that LDA is the first dimension of the actual array declaration, not the dimensions of the submatrix being operated on.

A secondary BLAS interface bug: the C CBLAS API and the Fortran BLAS API have different argument orders and character conventions. The Fortran DGEMM call takes the transpose flags as character arguments: TRANSA='N' for no transpose, TRANSA='T' for transpose. The exact character and case may be implementation-dependent in some BLAS libraries. A developer who calls the Fortran BLAS from C code via the CBLAS interface and uses the wrong row/column major flag will produce a transposed result. Retainer work that mixes Fortran and C numerical code needs to verify that the matrix layout convention is consistent at every BLAS or LAPACK call site.

Typical Fortran retainer work and what it looks like in a work log

Array bounds mismatch is the largest category of Fortran retainer work that produces wrong simulation output with no runtime error under production compiler settings. A subroutine that loops from 1 to SIZE on an assumed-shape argument that has a 0-based lower bound reads the wrong elements from the array. The results are wrong by one index in one dimension. For a climate simulation, this means the vertical temperature profile is shifted by one pressure level — the kind of error that produces scientifically plausible-looking output that passes casual inspection and fails only when compared against a reference solution. Work log entry: “compute_radiation: assumed-shape dummy temperature(:,:,:) inherited lower bound 0 from calling unit declaration temperature(0:nz-1,ny,nx); DO loop in subroutine started at index 1, skipping surface level 0; surface temperature not processed; vertical profile shifted 1 level; fix: changed main-program declaration to temperature(nz,ny,nx) default 1-based; reconciled all access sites; wrong simulation outputs before: 1; after: 0; 6h.”

OpenMP race condition is the second category. A loop accumulator declared SHARED instead of REDUCTION produces non-deterministic totals. Work log entry: “integrate_flux: global_total declared SHARED in !$OMP PARALLEL DO; 16 threads updated global_total concurrently without synchronization; total wrong and non-deterministic (different value each run); fix: added REDUCTION(+:global_total) to !$OMP PARALLEL DO directive; removed explicit global_total = global_total + thread_contribution line (handled by REDUCTION clause); wrong totals before: every multi-thread run; after: 0; 5h.”

BLAS leading dimension error is the third category. A DGEMM call with wrong LDA produces wrong matrix multiply results. Work log entry: “DGEMM call in solve_system: LDA passed as M=300 (submatrix row count); actual array A declared as A(1000, 800); BLAS read column offsets using LDA=300, not 1000; columns 2 and beyond started 700 rows early in memory; matrix multiply result numerically wrong for all off-diagonal elements; no BLAS error code; fix: changed LDA=M to LDA=1000 (declared array row dimension); wrong solve results before: all runs; after: 0; 4h.”

How HourTab tracks Fortran developer retainer hours

Fortran retainer work is invisible by the same mechanism that makes Fortran array operations efficient: the language does not add bounds checks by default in production builds. An assumed-shape argument that loops from 1 to SIZE on a 0-indexed array reads valid memory — the next element past the intended start is a real array element, not a segfault. The output is numerically plausible. The wrong value is one level off in a vertical dimension of a climate model — a difference that matches the noise level of the physical parameterization being maintained. The connection between a missing zero lower bound reconciliation and a shifted vertical profile requires knowing that assumed-shape dummy arguments inherit calling-unit bounds in Fortran 90 and later.

The work log needs to name the mechanism: the array declaration in the calling unit, the dummy argument type in the subroutine, which indexing convention each unit used, and how the discrepancy produced the wrong output. A log entry that says “fixed vertical profile bug, 6h” does not explain why the fix was a declaration change from REAL :: temperature(0:nz-1, ny, nx) to REAL :: temperature(nz, ny, nx). A log entry that names the assumed-shape inheritance mechanism and the loop start offset is auditable.

HourTab gives Fortran developers a public retainer-hours URL they send to clients — national laboratories maintaining climate models, university research groups maintaining fluid dynamics and quantum chemistry codes, defense contractors maintaining radar and signal processing simulations, and engineering firms maintaining finite element analysis programs. For Fortran retainers, each work log entry should name the Fortran mechanism: which array declaration, which dummy argument type, which loop bounds, which OpenMP clause, which BLAS parameter. Comparative context: Fortran retainer work has conceptual overlap with other numerical computing environments — COBOL (where PIC clause decimal alignment produces the analogous silent-precision-loss category of hard-to-log diagnostic work); Ada (where range-constrained subtypes provide the analogous explicit bounds that IMPLICIT NONE and assumed-shape bounds in Fortran approximate); and C (where 0-based array indexing is the default, creating the inverse off-by-one risk when C-background developers write Fortran). Fortran is uniquely positioned as the language for the numerical simulation codes that have been validated against physical measurements over decades and cannot be replaced without re-validation.

FAQ: Fortran developer retainers