Case Study 1: Porting a Legacy Trapezoid Integrator

"The old code was thirty lines and it was right. My job was to keep the second part true while fixing the first."

Executive Summary

The dictionary is only useful if you can apply it to code you did not write. This case study takes a small, complete FORTRAN 77 program — a trapezoidal-rule numerical integrator, the kind that appears in ten thousand engineering codes — and ports it to modern Fortran one construct at a time, using Chapter 19's entries as a checklist. The program is deliberately compact but it is not a toy: in thirty lines it packs a COMMON block, a DATA statement, a statement function, a GOTO accumulation loop, implicit typing, and fixed-form source — six dictionary entries in one place. We translate each, and then we do the thing that separates a port from a rewrite: we confirm, by hand, that the modern program computes exactly the same number as the original.

Skills applied: COMMON → module or arguments (§19.1); DATA → initializer (§19.1); statement function → internal pure function (§19.3); GOTO loop → do/exit (§19.2); implicit typing → implicit none + declarations (§19.4); fixed-form → free-form (§19.4); and the numerical-equivalence check from Chapter 18.

Background

You have inherited TRAP, which approximates $\int_a^b f(x)\,dx$ by the trapezoidal rule,

$$ \int_a^b f(x)\,dx \;\approx\; h\left(\tfrac{1}{2}f(a) + \sum_{k=1}^{n-1} f(a + kh) + \tfrac{1}{2}f(b)\right), \qquad h = \frac{b-a}{n}, $$

for the fixed integrand $f(x) = x^2$ on $[0, 3]$ with $n = 3$ intervals. Here is the code exactly as you found it — read it once through before we take it apart:

C     TRAP -- TRAPEZOIDAL-RULE INTEGRAL OF F(X) = X*X ON [A, B].
      PROGRAM TRAP
      COMMON /LIMITS/ A, B
      INTEGER N, K
      F(X) = X * X
      DATA A, B, N /0.0, 3.0, 3/
      H = (B - A) / REAL(N)
      AREA = 0.5 * (F(A) + F(B))
      K = 1
   10 CONTINUE
      X = A + REAL(K) * H
      AREA = AREA + F(X)
      K = K + 1
      IF (K .LE. N-1) GO TO 10
      AREA = AREA * H
      WRITE (*, 900) AREA
  900 FORMAT ('INTEGRAL = ', F6.2)
      END

Phase 1 — Inventory the Constructs

Before translating anything, list every legacy construct and the dictionary entry that handles it. This "parts list" is the whole method: you are not improvising, you are looking things up.

Note the one judgment call already visible in the table: COMMON /LIMITS/ can become a module, but a and b are nothing more than the integrator's inputs. The faithful-and-better move is to promote them to arguments of the integration routine, which is exactly the "how far past the literal translation to go" decision §19.1 flagged. A dictionary gives you the literal replacement; judgment tells you when a cleaner one is right there.

Phase 2 — Translate the Data and the Helper

Start with storage and the statement function, because they shape the interface. The integrand becomes an internal pure function with declared types (§19.3); the limits and interval count become arguments, not globals (§19.1); and the whole thing lives in a module so it has an explicit interface (§19.4, via Chapter 8):

module quad
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  private
  public :: trapezoid
contains

  pure function trapezoid(a, b, n) result(area)
    real(dp), intent(in) :: a, b        ! were COMMON /LIMITS/ -> now arguments
    integer,  intent(in) :: n
    real(dp) :: area, h, x
    integer  :: k
    h = (b - a) / real(n, dp)
    area = 0.5_dp * (f(a) + f(b))
    do k = 1, n - 1                      ! was the GOTO accumulation loop
       x = a + real(k, dp) * h
       area = area + f(x)
    end do
    area = area * h
  end function trapezoid

  pure function f(x) result(y)           ! was the statement function F(X) = X*X
    real(dp), intent(in) :: x
    real(dp) :: y
    y = x * x
  end function f

end module quad

Every legacy construct from the parts list is now handled: the COMMON globals are intent(in) arguments, the statement function is a pure module function, the GOTO loop is a counted do, and implicit none plus real(dp) declarations replace implicit typing. The pure attribute is honest here — the integrator has no side effects — and it lets the compiler inline f into the loop.

Phase 3 — Translate the Driver

The main program becomes a three-line driver that calls the routine and prints the result with an inline format (no separate numbered FORMAT statement needed):

program integrate
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use quad, only: trapezoid
  implicit none
  print '(a, f6.2)', 'integral = ', trapezoid(0.0_dp, 3.0_dp, 3)
end program integrate

$ gfortran -std=f2018 -Wall -O2 quad.f90 integrate.f90 -o integrate && ./integrate
integral =   9.50

Phase 4 — Prove the Numerics Did Not Change

This is the step that makes it a port and not a hopeful rewrite. Trace both programs by hand and confirm they produce the same number. With $a = 0$, $b = 3$, $n = 3$, the interval width is $h = (3-0)/3 = 1$.

• The endpoint contribution: $0.5\,(f(0) + f(3)) = 0.5\,(0 + 9) = 4.5$.
• The interior sum, $k = 1$ then $k = 2$: at $x = 1$, add $f(1) = 1$; at $x = 2$, add $f(2) = 4$. Running total $4.5 + 1 + 4 = 9.5$.
• Multiply by $h$: $9.5 \times 1 = 9.5$.

Both the legacy GOTO loop (which runs $k = 1$ while $k \le n-1 = 2$) and the modern do k = 1, n-1 execute for exactly $k = 1, 2$, so they accumulate the identical two terms. The result is 9.50 from both — a faithful translation. (For the record, the exact integral is $\int_0^3 x^2\,dx = 9$; the trapezoidal rule overestimates by $0.5$ because $x^2$ is convex. That error belongs to the method, not to the translation, and it is identical in both versions — which is the point.)

The check that matters: a modernization is only correct if the numbers survive it. Here they do, exactly, because every construct we changed was a matter of form — how the loop and the helper are written — not of arithmetic. When a translation also touches arithmetic (reordering a sum, changing precision), you drop from "bit-for-bit identical" to "equal within tolerance," the distinction Chapter 18 draws.

Phase 5 — What the Port Bought

Tally the improvements, because they justify the effort to whoever signs off on the change:

• The interface is explicit and checked. trapezoid(a, b, n) states its inputs; the compiler verifies every call. The old COMMON /LIMITS/ was invisible plumbing any routine could reach into.
• The integrand is reusable and testable. f is a real pure function you can call, test, and swap. To integrate a different function you now pass a different routine (an easy next step — see the Extensions); the statement function was welded in place.
• The loop is visible and safe. do k = 1, n-1 cannot fall through to the wrong label, and its bounds are in plain sight. implicit none guarantees a mistyped aera would be caught, not silently summed.
• The numbers are unchanged. The science the code encodes — the trapezoidal approximation — is byte-for-byte what it was. We improved the engineering and preserved the result, which is the entire ethic of Part IV.

Discussion Questions

1. We promoted COMMON /LIMITS/ A, B to arguments rather than to a module. Give one situation where a module would be the better choice, and one where arguments clearly win. What distinguishes them?
2. The legacy GOTO loop tests K .LE. N-1 at the bottom, so it always runs at least once. The modern do k = 1, n-1 runs zero times when n = 1. Does that difference change any result for $n \ge 2$? What about the edge case $n = 1$, and which behavior is correct for the trapezoidal rule?
3. f is marked pure. What specifically does that promise let the compiler do inside trapezoid, and why was the promise safe to make here but not for a version of f that printed a debug line?

Your Turn: Extensions

• Option A. Make the integrand a procedure argument: change trapezoid(a, b, n) to trapezoid(func, a, b, n) where func is passed in (an interface block or an abstract interface), and integrate both $x^2$ and $\sin x$ with the same routine. This is the modern answer to "the statement function was welded in."
• Option B. Add an intent(out) optional argument n_evals that reports how many times f was called, and confirm it equals $n + 1$. A textbook use of an optional output (Chapter 6).
• Option C. Replace the explicit accumulation do loop with a single whole-array expression: build the vector of sample points with an implied-do, apply f elementally, and sum with the sum intrinsic and a half-weight on the ends. Compare readability, and confirm the result is still 9.50.

Key Takeaways

• Inventory first, translate second. List every legacy construct and its dictionary entry before you change a line; the port becomes a lookup, not an improvisation.
• The literal replacement is a floor, not a ceiling. COMMON → module is faithful, but COMMON → arguments was better here. The dictionary tells you what is possible; judgment picks what is right.
• A port is defined by the numbers surviving it. Trace both versions by hand; if the result differs when you only changed form, you made a mistake, not a modernization.
• Six dictionary entries fit in thirty lines. Real legacy code is dense with translatable constructs, and the same handful of moves — COMMON, DATA, statement function, GOTO, implicit typing, fixed-form — covers the overwhelming majority of what you will meet.