Chapter 18: Modernizing Legacy Fortran — A Practical Migration Guide

"To me, legacy code is simply code without tests." — Michael Feathers, Working Effectively with Legacy Code

Overview

In Chapter 17 you learned to read old Fortran without panicking — to see past the fixed columns, the COMMON blocks, the GOTOs, and the variables typed by the first letter of their names, to the validated science underneath. This chapter is where you pick up the tools and change it. Not rewrite it — change it, one careful, reversible step at a time, until the crusty FORTRAN 77 relaxation kernel you met last chapter is clean, modern, modular Fortran that a newcomer can read and the compiler can check — and, crucially, computes exactly the same answer it always did.

That last clause is the whole discipline. Modernization is not an aesthetic exercise and it is not a rewrite; it is a sequence of behavior-preserving transformations, each small enough to verify, applied to a program that keeps working the entire time. The value in a legacy scientific code is not its source text — that part is genuinely unpleasant — but its validated numerics, the years of agreement with experiment that no rewrite can inherit. Our prime directive, the one this whole part of the book turns on, is therefore simple: improve the engineering without touching the science. By the end of the chapter you will have done exactly that to a real program, and you will be able to prove you did it.

In this chapter, you will learn to:

• Follow an eight-step recipe that takes any FORTRAN 77 program to modern Fortran in a defined order, from the safest, most mechanical changes to the most invasive.
• Work incrementally — keep the code compiling and passing its tests after every edit, so a mistake is caught in the step that caused it, not ten steps later.
• Build a regression test around a legacy program and use it to establish numerical equivalence, choosing correctly between demanding bit-for-bit agreement and accepting close enough.
• Map each legacy construct — COMMON, GOTO, EQUIVALENCE, statement functions, implicit typing, assumed-size arrays — to its modern replacement, and know why the replacement is better.
• Carry out a complete, worked modernization of the PLATE kernel, watching it transform before and after, and verifying the modern version reproduces the original's steady state.

Learning Paths

How to read this chapter by track. - 🔧 Legacy ("I inherited old code") — this is your chapter; read every section closely. §18.2 and §18.3 are the professional core: how to change code safely and how to prove you didn't break it. - 🔬 Scientist ("my old solver is slow") — you cannot profile or parallelize what you cannot build cleanly. Read §18.1 for the moves and §18.4 for the worked example; §18.3 is how you protect your validated results while you improve the code. - 📖 Standard — read straight through; the before/after in §18.4 shows the modern features of Parts I–II doing real work against their legacy counterparts. - ⚡ HPC — modernization is the prerequisite to parallelization. A clean modular code with intent and no aliasing is what OpenMP and MPI (Part VIII) can actually get a grip on; skim §18.1, focus on §18.4.

18.1 The Eight-Step Recipe

Faced with a few thousand lines of FORTRAN 77, the temptation is to start typing and fix whatever offends you first. Resist it. Modernization goes better as an ordered sequence, because the early steps make the later ones safe. implicit none (step 1) exposes the variables you will need to understand before you can split code into modules (step 3); modules give you the explicit interfaces that make adding intent (step 4) meaningful; and so on. The order below runs, roughly, from the safest and most mechanical changes to the most structural. It is a default, not a law — on a real code you will interleave and prioritize by payoff — but when in doubt, go in this order.

Throughout, we lean on one callout more than any other:

🔧 Modern vs Legacy: the whole chapter is a sustained use of this device. For each step you will see the FORTRAN 77 on the left, in spirit, and the modern Fortran on the right — and a sentence on why the modern form is not merely prettier but safer, faster to reason about, or checkable by the compiler.

Before we start, name the activity honestly. Each step is a refactoring: a change to the form of the code that leaves its observable behavior unchanged. The word matters because it sets the acceptance test. A refactoring that changes the answer is not a refactoring — it is a bug, and §18.3 is how you catch it.

Step 1 — Add implicit none

The single most valuable line you can add to old Fortran is the one that turns off implicit typing (the convention, defined in Chapter 17, that undeclared names starting with I–N are integers and everything else is real). With implicit none in force, every variable must be declared, and a mistyped name becomes a compile error instead of a silently created new variable holding garbage.

🔧 Modern vs Legacy: implicit typing hides one of the most expensive bugs in scientific computing.

fortran C Legacy: no declarations. A typo creates a new zero-ish variable, silently. VELOCITY = 0.0 KINETIC = 0.5 * MASS * VELOCTY**2 C VELOCTY (typo) is implicitly REAL, undefined -> wrong energy, no warning.

fortran ! Modern: implicit none makes the typo a compile-time error. real(dp) :: velocity, kinetic, mass velocity = 0.0_dp kinetic = 0.5_dp * mass * velocty**2 ! error: 'velocty' has no IMPLICIT type

Adding implicit none to a unit that relied on implicit typing will not compile until you have declared every variable — which is the point. The compiler hands you a checklist of every name in the routine, and working through it is how you learn the routine. This is the step that turns reading into understanding.

Step 2 — Convert fixed-form to free-form

FORTRAN 77 is written in fixed-form source (Chapter 17): columns 1–5 for statement labels, column 6 for continuation, code confined to columns 7–72, a C in column 1 marking a comment. Modern Fortran is free-form — statements go anywhere, ! starts a comment, & continues a line, and the compiler stops caring what column you are in. The conversion is almost entirely mechanical, and you should let a tool do the bulk of it.

💡 Intuition: think of this as changing the file format, not the program. Nothing about the logic changes; you are lifting the code out of the punch-card grid it was born in.

Tools such as findent and fprettify convert and reindent fixed-form to free-form automatically; the open-source plusFORT and academic CamFort go further. Run the tool, then read the diff — automated conversions are reliable but not infallible (continuation lines and long statements are where they occasionally slip), and reviewing the change is itself part of learning the code.

📜 From History: the 72-column limit is the width of an 80-column punch card minus the 8 columns reserved for a card sequence number, so that a dropped deck could be re-sorted. Your code has been shaped for decades by a physical object most working programmers have never held.

Step 3 — Replace COMMON blocks with modules

This is the central move of modernization, and the one Chapter 8 spent a whole chapter preparing you for. A COMMON block is global state by memory overlay: a set of variables laid end to end in a shared region, with every program unit that names the block trusting — without any check — that it lists the same variables, in the same order, with the same types. Get the list wrong in one routine and you reinterpret the bytes, silently and catastrophically.

A module replaces that with a real namespace: named, typed variables that every user uses by name, with the compiler enforcing agreement.

🔧 Modern vs Legacy: the same shared grid, two eras.

fortran C Legacy: COMMON overlays memory; every unit must repeat this line EXACTLY. PARAMETER (NMAX = 21) COMMON /GRID/ T(NMAX,NMAX), N

fortran ! Modern: a module. Declared once; used by name; checked by the compiler. module grid_data implicit none integer, parameter :: dp = selected_real_kind(15, 307) integer, parameter :: nmax = 21 real(dp) :: t(nmax, nmax) integer :: n end module grid_data

The mechanical translation — one module per COMMON block, one module variable per block member — is a faithful, behavior-preserving first move: it keeps the global nature of the data (which the rest of the code still depends on) while making it typed and checked. Removing the globality altogether, by passing the data as arguments, is a later and larger refinement (steps 4–5). Do the safe thing first.

🚪 Threshold Concept — modules replace the global mess of COMMON. Once the shared state of a program lives in modules rather than COMMON blocks, the compiler can see the whole data-flow: who declares what, who reads it, who writes it. Bugs that COMMON could only produce at run time (a mismatched block, a type pun) become impossible or compile-time errors. This one change does more for the reliability of an old code than any other, which is why it sits at the center of the recipe.

Step 4 — Add intent to every dummy argument

Once routines pass data as arguments rather than sharing it through COMMON, declare the direction of every argument with intent: intent(in) for data the routine only reads, intent(out) for data it produces, intent(inout) for data it modifies in place. FORTRAN 77 had no such concept — every argument was silently inout, and passing a constant to a routine that overwrote it was a memorable way to change the value of the literal 2 for the rest of your program.

🔧 Modern vs Legacy: intent is a contract the compiler enforces.

fortran C Legacy: no way to say "SUBROUTINE only reads A". Everything is modifiable. SUBROUTINE SCALE(A, N, S) DIMENSION A(N)

fortran ! Modern: the interface documents and enforces who may change what. subroutine scale(a, s) real(dp), intent(inout) :: a(:) ! scaled in place real(dp), intent(in) :: s ! read only; compiler forbids writing to s

intent is documentation the compiler checks, an optimization hint, and a bug-catcher in one. It is also free: adding it costs a few keystrokes per argument and can only make the code safer.

Step 5 — Replace assumed-size arrays with assumed-shape

Legacy code passes arrays as assumed-size — real a(*) or real a(n) with n passed alongside — which tells the callee nothing about the array's real bounds; the routine trusts the caller's separately-passed dimensions and walks off the end if they disagree. Modern code uses assumed-shape arrays, real(dp), intent(in) :: a(:,:), where the array carries its own shape and the callee asks it with size and shape (Chapter 6).

🔧 Modern vs Legacy: the array should know its own size.

fortran C Legacy: dimensions travel separately; a mismatch is undetected memory corruption. SUBROUTINE RESID(A, NX, NY, R) DIMENSION A(NX, NY)

fortran ! Modern: shape travels with the array; query it where you need it. subroutine resid(a, r) real(dp), intent(in) :: a(:,:) real(dp), intent(out) :: r integer :: nx, ny nx = size(a, 1); ny = size(a, 2)

Assumed-shape also enables bounds checking (-fcheck=all) to catch an out-of-range index at run time — a safety net the separately-passed-dimension style makes impossible.

Step 6 — Replace GOTO with structured control flow

The GOTO, computed GOTO, and arithmetic IF of FORTRAN 77 (Chapter 17) become the structured constructs of Chapter 4: do loops with exit and cycle, if … then … else, and select case. The classic pattern — a GOTO-based convergence loop — maps directly onto a do loop with an exit:

🔧 Modern vs Legacy: the same iterate-until-converged loop, unknotted.

fortran C Legacy: a backward GOTO forms the loop; the exit test is a forward GOTO. ITER = 0 10 CONTINUE ITER = ITER + 1 C ... do one sweep, compute DMAX ... IF (DMAX .GT. TOL .AND. ITER .LT. MAXIT) GO TO 10

fortran ! Modern: the loop structure is visible; no labels, no jumps. iters = 0 do iters = iters + 1 ! ... do one sweep, compute dmax ... if (dmax <= tol .or. iters >= maxit) exit end do

Note the exit condition is the logical negation of the legacy continue condition — "keep going while DMAX > TOL and iterations remain" becomes "stop when dmax <= tol or we hit the cap." Getting that negation right is exactly the sort of thing a regression test (§18.3) exists to confirm.

Step 7 — Replace EQUIVALENCE with types or transfer

EQUIVALENCE forces two names to share storage — the most dangerous idiom in FORTRAN 77, used variously to save memory (overlay a scratch array), to reinterpret bits (view a real as an integer), or to alias a 2-D array as 1-D. Each use has a clean modern replacement, and none of them aliases memory behind the compiler's back.

🔧 Modern vs Legacy: replace memory tricks with intent-revealing operations.

fortran C Legacy: alias the 2-D grid as a 1-D vector to take a whole-array max. DIMENSION T(NMAX,NMAX), TVEC(NMAX*NMAX) EQUIVALENCE (T(1,1), TVEC(1)) TMAX = TVEC(1) DO 10 K = 2, N*N IF (TVEC(K) .GT. TMAX) TMAX = TVEC(K) 10 CONTINUE

fortran ! Modern: no aliasing needed -- ask the array directly. tmax = maxval(t(1:n, 1:n))

When the aliasing was a genuine reinterpretation of bits — the one use with no array-intrinsic replacement — reach for the transfer intrinsic, which converts a value's bit pattern to another type explicitly and visibly, so the reinterpretation is a documented operation rather than a hidden overlay. When EQUIVALENCE bundled related quantities into one block of storage, a derived type expresses the grouping with names and types instead.

Step 8 — Add error handling

FORTRAN 77 codes classically fail by marching off the end of an array, dividing by a zero nobody checked for, or reading a malformed input file and continuing with rubbish. Modern Fortran gives you the tools to fail loudly and early (Chapter 13): error stop with a message for an unrecoverable condition, stat= and errmsg= on allocate, iostat= and iomsg= on I/O, and explicit validation of inputs at the point of entry.

🔧 Modern vs Legacy: validate, then compute.

fortran C Legacy: trust the input; a bad N corrupts memory downstream, silently. READ (5,*) N

fortran ! Modern: check the precondition; stop with a diagnosis if it fails. read (unit, *, iostat=ios) n if (ios /= 0) error stop 'plate: could not read grid size' if (n < 3 .or. n > nmax) error stop 'plate: grid size out of range'

This step is last not because it is least important but because it is easiest to add once the code is clean: with modules, intent, and assumed-shape in place, you can see where the inputs enter and where the assumptions live, and put a guard on each.

🔄 Check Your Understanding. (1) Why does implicit none come first in the recipe? (2) Which single step does the most for reliability, and why? (3) In step 6, the modern exit condition is the negation of the legacy GOTO condition — what property of the migration must you verify because of that?

Answers

(1) It forces you to declare — and therefore understand — every variable, which every later step depends on; it also immediately catches typo-bugs. (2) Replacing COMMON with modules: it turns unchecked global memory overlay into compiler-checked, typed, named state, eliminating a whole class of run-time-only bugs. (3) That the negation is exact — that the loop stops on precisely the same iteration it always did — which is what a regression test confirms.

18.2 Modernizing Incrementally Without Breaking a Working Code

The recipe tells you what to change. This section is about how to change it without ever holding a broken program in your hands. The method has a name.

Definition (incremental modernization). Transforming a working program through a series of small, individually verified changes — compiling and re-running its tests after each one — so that the code remains correct and buildable at every step, and any regression is localized to the single edit that caused it. It is the opposite of a "big-bang" rewrite, in which the program does not run again until a large body of new code is finished and all its bugs arrive at once.

The case for incrementalism is a case about where bugs come from and how you find them. If you make one behavior-preserving change and the regression test still passes, you have strong evidence that change was safe. If you make forty changes and the test fails, you have forty suspects and no alibi for any of them. The cost of finding a bug scales horribly with the number of unverified changes between you and it — so you drive that number to one.

🚪 Threshold Concept — never rewrite what you can refactor. The instinct, facing ugly code, is to start over with a blank file. For validated scientific software this instinct is almost always wrong, and understanding why changes how you work forever. A rewrite throws away the one irreplaceable asset — the validated behavior — and does not get it back until the rewrite is complete, tested, and re-validated, which for real codes routinely takes longer than the original did and reproduces only a fraction of its capability. Refactoring keeps the validated behavior the whole time and improves the engineering around it. The working program is never out of your hands. Tape the motto to your monitor: never rewrite what you can refactor.

A concrete loop for a single increment looks like this:

  1. Pick the smallest useful change (one routine's implicit none; one COMMON block).
  2. Make it.
  3. Compile.  Does it build clean, with -Wall?
  4. Run the regression test.  Does the output still match the reference?
  5. Commit, with a message naming the step.  (git -- Chapter 37.)
  6. Repeat.

Two supports make this loop safe. The first is version control. Every green step is a commit; if a later change breaks the test, git diff against the last green commit shows you the exact lines to suspect, and git revert gets you back to working code instantly. Modernization without version control is mountaineering without a rope. (We treat git for scientific code properly in Chapter 37.)

The second support is the regression test itself, and there is a subtlety about ordering it with respect to the work. You must write the test before you change code you do not fully understand. A test written against the current, trusted behavior — even behavior you cannot yet explain — pins that behavior so you can refactor underneath it. This is Michael Feathers' idea of a characterization test: a test whose job is not to specify what the code should do but to capture what it currently does, so that any change you make is measured against the code's own established behavior. For a legacy scientific code, the "specification" often exists only as the running program; the characterization test makes that specification executable.

⚠️ Common Pitfall: "I'll modernize the whole file, then test." This is the big-bang rewrite wearing a refactoring costume, and it fails the same way: when the test finally goes red, the cause is buried somewhere in a hundred edits. Keep the granularity brutal — one conceptual change per test run. It feels slow. It is dramatically faster than debugging a heap of simultaneous changes.

Not every step is equally risky, and it helps to know which are which. Steps 1, 2, and 6 (implicit none, free-form, GOTO → do) are local and mechanical — confined to one routine, changing form not logic, and largely tool-assisted. Steps 3, 4, and 5 (COMMON → modules, intent, assumed-shape) are structural: they change interfaces, so they ripple across every caller, and they are where you go slowest and lean hardest on the test. Step 7 (EQUIVALENCE) is case by case and demands the most thought, because you must first work out which of its several purposes the original author intended. Sequence your work to get the cheap, safe wins first: they make the code readable, which makes the risky steps less risky.

🔗 Connection: the incremental method is not special to Fortran — it is the mainstream discipline of software refactoring, adapted to numerical code. What is special to numerical code is the acceptance test, because "the same behavior" means "the same numbers," and numbers have tolerances. That is the subject of the next section.

18.3 Testing the Migration: Regression Tests and Numerical Equivalence

You cannot refactor safely without a way to tell, after each change, whether the program still does what it did. For numerical code that means comparing its output numbers against a trusted reference.

Definition (regression test). A test that runs the program on a fixed input and compares its output against a stored reference — a "golden" output captured earlier from the trusted version — reporting failure if they differ by more than an allowed amount. Its purpose is to catch a regression: a change that alters behavior that was previously correct. In a migration, the reference is the output of the original, validated program, and the regression test is what lets you claim your modern version still computes the same result.

The reference is captured once, from the original code, before you touch anything:

$ gfortran -std=legacy plate-legacy.f -o plate-legacy
$ ./plate-legacy > reference.txt        # the golden output, from the trusted code

From then on, every modernization step is followed by rebuilding the modern version, running it, and comparing its output to reference.txt. The comparison is the delicate part, and it turns on a definition.

Definition (numerical equivalence). Two programs are numerically equivalent on a given input when they produce the same results — either bit-for-bit identical (every floating-point value equal to the last bit) or identical to within a stated tolerance (every value agreeing to some absolute or relative error you decide is acceptable). Which standard you can demand depends entirely on whether your changes altered the arithmetic.

The distinction is the heart of testing a migration, so make it sharp:

Bit-for-bit reproducibility is the strongest claim: the modern program's every output value equals the legacy program's to the last bit. You can demand it — and should — whenever your changes did not touch the arithmetic. The mechanical and structural steps of §18.1 are all arithmetic-preserving: implicit none, free-form conversion, COMMON → modules, intent, assumed-shape, GOTO → do. None of them changes which operations happen, in which order, at which precision, so none of them changes the rounding, so the results are identical to the bit. If a "pure refactoring" step breaks bit-for-bit agreement, you did not do a pure refactoring — you introduced a change, and the test just caught it. That is the test working exactly as intended.

"Close enough" (tolerance-based) equivalence is what you fall back to when you deliberately change the numerics — and some modernizations do. Promote a single-precision real to real(dp) (Chapter 3) and every value shifts in its low-order bits: that is an improvement, not a regression, but bit-for-bit is gone and you must validate to a tolerance instead. Reorder a summation, let the compiler contract a multiply-add into a fused fma, or turn on -ffast-math, and the same thing happens (Chapter 20 explains why floating-point arithmetic is not associative). When you change the math on purpose, choose a tolerance from the problem's accuracy needs and compare against that.

⚠️ Common Pitfall — bit-for-bit is a property of the build, not just the source. Even identical source can produce different bits under different compilers, optimization levels, or fast-math flags, because those change instruction selection and expression association. If you promise a collaborator "bit-for-bit reproducible," you are promising a specific compiler and a specific set of flags — record them (Chapter 30). Reproducibility that is not pinned to a build is a wish, not a guarantee.

There is one more trap, and it is the most common of all: comparing output as text.

⚠️ Common Pitfall: a naïve diff reference.txt modern.txt compares characters, and will flag a migration as "broken" over changes that are not numerical at all — a legacy WRITE (6,'(I5)') that padded an integer to five columns versus a modern print '(i0)' that does not, a D exponent versus an E, a trailing space. These are formatting differences, not numerical ones. A regression test for numerical code must parse the numbers out and compare them within tolerance, ignoring whitespace and format. Compare the physics, not the pretty-printing.

We will meet the industrial-strength version of all this — the pFUnit unit-testing framework, regression suites, and continuous integration across multiple compilers — in Chapter 37. Here, the tiny hand-rolled harness is enough to make the point and to guard the PLATE migration in §18.4: capture a reference, and after every step confirm the numbers still match.

🧩 Try It Yourself: before reading the worked example, predict which category the PLATE modernization will fall into. We keep double precision throughout and never reorder the four-neighbor average — so should we expect bit-for-bit agreement, or only close enough? Write down your answer; §18.4 confirms it (and shows that for the small test grid the agreement is not merely bit-for-bit but exact, with no rounding at all).

18.4 A Full Worked Modernization of PLATE

Now we do it for real. This is the program you met in Chapter 17: PLATE, a FORTRAN 77 code that finds the steady-state temperature of a square plate by Jacobi relaxation of Laplace's equation. Physically, the plate's edges are held at fixed temperatures and we want the interior temperature once it stops changing; numerically, each interior point is repeatedly replaced by the average of its four neighbors until the field settles. At steady state every interior value satisfies the discrete Laplace equation

$$ T_{i,j} = \tfrac{1}{4}\left(T_{i-1,j} + T_{i+1,j} + T_{i,j-1} + T_{i,j+1}\right), $$

which is the five-point stencil with no source term. (The physics and the stencil are developed properly in Chapter 24; here the program is a legacy artifact to modernize, not a numerical-methods lesson.)

The "before": legacy PLATE in full

Here is the code, fixed-form and unmodernized. Read it as archaeology — every legacy trait from Chapter 17 is on display: IMPLICIT DOUBLE PRECISION, a COMMON /GRID/ block repeated in every unit, a DATA statement for the edge temperatures, a statement function for the update, and a GOTO-based convergence loop.

C=====================================================================
C  PLATE -- steady-state temperature of a square plate.
C           2-D Laplace equation, solved by Jacobi relaxation.
C           FORTRAN 77, fixed-form.  *** LEGACY CODE -- do not imitate. ***
C           Compile:  gfortran -std=legacy plate-legacy.f -o plate-legacy
C=====================================================================
      PROGRAM PLATE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NMAX = 21)
      COMMON /GRID/ T(NMAX,NMAX), N
C     Grid size (a 4x4 test plate; interior is 2x2).
      N = 4
      CALL SETBC
      CALL RELAX
      CALL OUTPT
      STOP
      END
C
      SUBROUTINE SETBC
C     Boundary conditions: top edge hot, other three edges cold.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NMAX = 21)
      COMMON /GRID/ T(NMAX,NMAX), N
      DATA THOT, TCOLD / 100.0D0, 0.0D0 /
      DO 20 J = 1, N
        DO 10 I = 1, N
          T(I,J) = 0.0D0
   10   CONTINUE
   20 CONTINUE
      DO 30 I = 1, N
        T(I,1) = TCOLD
        T(I,N) = TCOLD
   30 CONTINUE
      DO 40 J = 1, N
        T(N,J) = TCOLD
   40 CONTINUE
      DO 50 J = 1, N
        T(1,J) = THOT
   50 CONTINUE
      RETURN
      END
C
      SUBROUTINE RELAX
C     Jacobi relaxation with a GOTO convergence loop.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NMAX = 21)
      COMMON /GRID/ T(NMAX,NMAX), N
      DOUBLE PRECISION TNEW(NMAX,NMAX)
      DATA TOL, MAXIT / 1.0D-6, 1000 /
C     Statement function: four-neighbour average.
      AVG(TL,TR,TB,TA) = 0.25D0 * (TL + TR + TB + TA)
      ITER = 0
   60 CONTINUE
      ITER = ITER + 1
      DMAX = 0.0D0
C     One Jacobi sweep: TNEW from the OLD T (interior only).
      DO 80 J = 2, N-1
        DO 70 I = 2, N-1
          TNEW(I,J) = AVG(T(I-1,J), T(I+1,J), T(I,J-1), T(I,J+1))
          DIFF = ABS(TNEW(I,J) - T(I,J))
          IF (DIFF .GT. DMAX) DMAX = DIFF
   70   CONTINUE
   80 CONTINUE
C     Copy the interior back into T.
      DO 100 J = 2, N-1
        DO 90 I = 2, N-1
          T(I,J) = TNEW(I,J)
   90   CONTINUE
  100 CONTINUE
      IF (DMAX .GT. TOL .AND. ITER .LT. MAXIT) GO TO 60
      WRITE (6,'(A,I5,A)') ' converged in ', ITER, ' iterations'
      RETURN
      END
C
      SUBROUTINE OUTPT
C     Print the interior temperatures.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NMAX = 21)
      COMMON /GRID/ T(NMAX,NMAX), N
      DO 120 I = 2, N-1
        WRITE (6,'(20F9.4)') (T(I,J), J = 2, N-1)
  120 CONTINUE
      RETURN
      END

Compile and run (this is our reference capture — the only time we build the legacy code):

console $ gfortran -std=legacy plate-legacy.f -o plate-legacy && ./plate-legacy > reference.txt $ cat reference.txt converged in 25 iterations 37.5000 37.5000 12.5000 12.5000

Take a moment with that output, because it is the fixed point of everything that follows. On this 4×4 grid — a top edge at 100°, three edges at 0° — the interior settles to two upper cells at 37.5° and two lower cells at 12.5°, reached in 25 Jacobi sweeps at a tolerance of $10^{-6}$. We will now transform the program step by step, and after every step this output must not move.

Verifying the reference by hand

We never take a program's word for its output; we compute it. With one hot edge and three cold, the two upper interior cells are equal by symmetry (call the value $a$) and the two lower cells are equal (call it $b$). Steady state means each equals the average of its four neighbors:

$$ a = \tfrac{1}{4}(100 + b + 0 + a) \;\Rightarrow\; 3a - b = 100, \qquad b = \tfrac{1}{4}(a + 0 + 0 + b) \;\Rightarrow\; a = 3b. $$

Substituting gives $8b = 100$, so $b = 12.5$ and $a = 3b = 37.5$ — exactly the printed field. The Jacobi sweeps march toward it from an all-zero interior:

  sweep   a = T(2,2)    b = T(3,2)     dmax (max change)
    1       25.0000       0.0000        25.0000
    2       31.2500       6.2500         6.2500
    3       34.3750       9.3750         3.1250
    4       35.9375      10.9375         1.5625
    5       36.7188      11.7188         0.7813
   ...          ...          ...          (dmax halves each sweep)
   25       37.5000      12.5000        ~7.5e-07  < 1.0e-06  -> stop

After the first sweep the maximum change halves every iteration, so it falls below $10^{-6}$ for the first time on sweep 25 — which is why the legacy code reports 25. Every value in this table is a dyadic fraction representable exactly in IEEE double precision, so no rounding occurs anywhere in this computation. That is worth flagging now: it means the modern version, doing the same operations in double precision, cannot merely round the same — it computes the identical values, exactly. This migration will be bit-for-bit, and on this grid, exact.

Applying the eight steps

Steps 1 and 2 (implicit none, free-form). Convert each unit to free-form and add implicit none, which forces us to declare every variable — the loop counters i, j, the scalars dmax, diff, tol, the integer iter, maxit, the arrays. Nothing about the logic changes; the program simply stops relying on the A–H, O–Z typing rule and states its types.

Step 3 (COMMON → module). The COMMON /GRID/ block — repeated verbatim in all four program units — collapses to a single grid_data module declared once and used where needed. This is the moment the program stops being a set of routines that happen to agree on a memory layout and becomes a set of routines that share a named, typed object.

Steps 4 and 5 (intent, assumed-shape). Here we make the larger design improvement the earlier steps set up: rather than communicating through global state at all, the relaxation becomes a subroutine that takes the temperature field as an assumed-shape intent(inout) argument. The global COMMON grid was faithful to the original, but passing the field explicitly is better still — the data-flow is now visible in the interface, the routine is reusable on any grid, and the array carries its own shape. We keep only the kind parameter dp in the module.

Step 6 (GOTO → do). The 60 CONTINUE … GO TO 60 loop becomes a clean do … exit … end do, with the exit condition the exact negation of the legacy continue condition, as discussed in §18.1.

Step 7 (EQUIVALENCE). PLATE uses none — a reminder that the recipe is a menu, not a fixed march; you apply the steps a given code actually needs. (You practice the EQUIVALENCE step on other code in the exercises.)

Step 8 (error handling). We guard the one precondition the arithmetic depends on: the grid must be square and at least 3×3 (or there is no interior to relax). A violation is an error stop with a message, not silent memory corruption.

The "after": modern plate.f90

! plate.f90 -- steady-state temperature of a square plate by Jacobi relaxation.
! The MODERNIZED form of plate-legacy.f -- free-form, implicit none, a module
! instead of COMMON, intent + assumed-shape arguments, a structured loop, error stop.
! Compile:  gfortran -std=f2018 -Wall -O2 plate.f90 -o plate && ./plate

module plate
  implicit none
  private
  integer, parameter, public :: dp = selected_real_kind(15, 307)   ! == real64
  public :: set_boundary, relax, print_interior
contains

  !> Dirichlet boundary conditions: top edge hot, the other three cold.
  subroutine set_boundary(t, thot, tcold)
    real(dp), intent(out) :: t(:,:)
    real(dp), intent(in)  :: thot, tcold
    integer :: n
    n = size(t, 1)
    t       = 0.0_dp       ! interior initial guess = 0
    t(:, 1) = tcold        ! left edge
    t(:, n) = tcold        ! right edge
    t(n, :) = tcold        ! bottom edge
    t(1, :) = thot         ! hot top edge (corners unused by the stencil)
  end subroutine set_boundary

  !> Jacobi relaxation to steady state; returns the sweep count in iters.
  subroutine relax(t, tol, maxit, iters)
    real(dp), intent(inout) :: t(:,:)
    real(dp), intent(in)    :: tol
    integer,  intent(in)    :: maxit
    integer,  intent(out)   :: iters
    real(dp) :: tnew(size(t,1), size(t,2)), dmax
    integer  :: i, j, n
    n = size(t, 1)
    if (n < 3 .or. size(t,2) /= n) error stop 'relax: grid must be square, n >= 3'
    tnew = t                                   ! keep boundary values fixed
    iters = 0
    do
      iters = iters + 1
      dmax = 0.0_dp
      do j = 2, n-1                            ! one Jacobi sweep from the OLD t
        do i = 2, n-1
          tnew(i,j) = 0.25_dp * (t(i-1,j) + t(i+1,j) + t(i,j-1) + t(i,j+1))
          dmax = max(dmax, abs(tnew(i,j) - t(i,j)))
        end do
      end do
      t = tnew
      if (dmax <= tol .or. iters >= maxit) exit
    end do
  end subroutine relax

  !> Print the interior temperatures, one grid row per line.
  subroutine print_interior(t)
    real(dp), intent(in) :: t(:,:)
    integer :: i, j, n
    n = size(t, 1)
    do i = 2, n-1
      write(*, '(*(f9.4))') (t(i,j), j = 2, n-1)
    end do
  end subroutine print_interior

end module plate

program plate_modern
  use plate, only: dp, set_boundary, relax, print_interior
  implicit none
  integer, parameter :: n = 4
  real(dp) :: t(n, n)
  integer  :: iters
  call set_boundary(t, thot = 100.0_dp, tcold = 0.0_dp)
  call relax(t, tol = 1.0e-6_dp, maxit = 1000, iters = iters)
  print '(a, i0, a)', 'converged in ', iters, ' iterations'
  call print_interior(t)
end program plate_modern

$ gfortran -std=f2018 -Wall -O2 plate.f90 -o plate && ./plate
converged in 25 iterations
  37.5000  37.5000
  12.5000  12.5000

Compare it to the reference. The interior field — 37.5000 37.5000 / 12.5000 12.5000 — is identical to the last digit, and (as we proved by hand) identical to the last bit: the migration is numerically equivalent, exactly. The program converges in the same 25 sweeps because the update expression, the operand order, the precision, the convergence measure, and the stopping rule were all preserved through every step.

🔧 Modern vs Legacy — read the two programs side by side. The modern version is shorter, and every line is checked in ways the original could not be: implicit none guarantees no undeclared variable silently exists; the module means no two routines can disagree about the grid; intent(in) on the boundary temperatures means relax physically cannot corrupt them; assumed-shape means the field carries its own size; error stop means a bad grid dies with a diagnosis instead of walking off the array. And it runs the identical arithmetic. That is the promise of the whole part made concrete: the science preserved, the engineering transformed.

⚠️ Common Pitfall — the log line is not the physics. Notice the one visible difference: the legacy WRITE (6,'(A,I5,A)') prints converged in 25 iterations (the I5 pads to five columns) while the modern print '(a,i0,a)' prints converged in 25 iterations (no padding). A character-level diff would flag the migration as changed; a numerical regression test parses the integer 25 from each and sees they match. This is exactly the §18.3 warning in the flesh: compare the numbers, not the whitespace.

🔄 Check Your Understanding. (1) Why is the PLATE migration bit-for-bit rather than only "close enough"? (2) What would break bit-for-bit agreement — name two changes a well-meaning modernizer might make. (3) The modern relax takes the field as an argument, but §18.1 step 3 put it in a module. Are both "correct" modernizations?

Answers

(1) We preserved precision (double), the update expression, the operand order, and the stopping rule, so no rounding changed. (2) Promoting the reals to a different precision; reordering or reassociating the four-neighbor sum; enabling -ffast-math/FMA contraction; switching the loop to Gauss–Seidel by updating t in place. Any of these changes the arithmetic. (3) Yes — the module-global form is the faithful mechanical step-3 result; passing the field as an intent(inout) assumed-shape argument is the further refinement of steps 4–5. Both preserve the answer; the argument form is the better final design.

Project Checkpoint

The running heat solver's Part IV side quest is the modernization of this very kernel, and your checkpoint is to take the first two steps yourself — the safe, high-value ones — and prove they changed nothing.

Your task: starting from plate-legacy.f, apply step 1 (implicit none, free-form) and step 3 (COMMON /GRID/ → a grid_data module) to produce an intermediate plate-ckpt.f90 that still uses the legacy GOTO loop and module-global data — steps 4–8 come later. Capture the legacy reference first, then confirm the intermediate reproduces it. The module and the driver's shape:

module grid_data
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: nmax = 21
  real(dp) :: t(nmax, nmax)     ! was:  COMMON /GRID/ T(NMAX,NMAX), N
  integer  :: n
end module grid_data

program plate_ckpt
  use grid_data
  implicit none
  real(dp) :: tnew(nmax, nmax), tol, dmax, diff
  integer  :: i, j, iter, maxit
  ! ... set boundary, then the SAME Jacobi loop, still with GO TO 100 ...
end program plate_ckpt

The full intermediate is in code/project-checkpoint.f90. It must still print, to the last digit:

$ gfortran -std=f2018 -Wall project-checkpoint.f90 -o plate-ckpt && ./plate-ckpt
converged in 25 iterations
  37.5000  37.5000
  12.5000  12.5000

Why this matters for the capstone. The Chapter 38 solver you are building is born modern, but real scientific software is not — it is inherited, and the skill of taking a validated legacy kernel to a clean, testable, parallelizable modern form is exactly what turns a solver into a maintainable code. You have now done, on a small kernel, what Chapter 37 will formalize as regression testing: change the form, hold the numbers. In Chapter 19 you finish the job, using the translation dictionary to carry the kernel through the remaining steps.

Summary

Modernizing legacy Fortran is a disciplined refactoring — improve the engineering, preserve the science, prove it — not a rewrite.

The two things worth memorizing: first, the recipe's order — implicit none, free-form, modules, intent, assumed-shape, structured control, kill EQUIVALENCE, add error handling. Second, the acceptance test that makes it safe — after every step, the numbers must match the reference, bit-for-bit when you did not touch the arithmetic.

Spaced Review

Retrieval practice on the chapters this one builds directly on: Chapter 8 (modules, the target of step 3) and Chapter 17 (the legacy constructs we are replacing).

1. (Ch. 8) Modernization step 3 replaces a COMMON block with a module. Name two things the compiler can check about a module's shared variables that it cannot check about a COMMON block's.
   

   Answer

   That every user refers to the variables by the same name with the same type and shape (a `use` cannot silently disagree the way a re-listed `COMMON` can), and that access respects `public`/`private`. The compiler also generates and checks explicit interfaces for module procedures, which `COMMON`-and-external-subroutine code lacks.

2. (Ch. 8) After step 3 you have module-global data; after steps 4–5 the data is passed as arguments. What does a module still legitimately own in the final plate.f90?
   

   Answer

   The kind parameter `dp` and the procedures themselves (`set_boundary`, `relax`, `print_interior`) — a module as a namespace of related procedures and constants, with the mutable field passed explicitly rather than shared globally.

3. (Ch. 17) In the legacy RELAX, the line AVG(TL,TR,TB,TA) = 0.25D0*(TL+TR+TB+TA) is a statement function. What is a statement function, and what modern construct replaces it?
   

   Answer

   A one-line function defined among the specification statements of a program unit, usable only within that unit. Modern Fortran replaces it with an internal procedure (a `function` inside `contains`), or — as we did — by inlining the tiny expression, since it was used once.

4. (Ch. 17) The legacy code writes IMPLICIT DOUBLE PRECISION (A-H,O-Z). Which variables in RELAX are integers under this rule, and how do you know?
   

   Answer

   `I`, `J`, `ITER`, and `MAXIT` — names beginning with a letter in the range `I`–`N`, which the implicit rule leaves as `INTEGER`; the `(A-H,O-Z)` clause makes only the *other* letters `DOUBLE PRECISION`. `T`, `TNEW`, `TOL`, `DMAX`, `DIFF`, `AVG` all begin outside `I`–`N` and are double.

5. (Ch. 8 + 17) Why does converting the four repeated COMMON /GRID/ declarations into one module reduce the chance of a bug, beyond being less typing?
   

   Answer

   Because four hand-copied `COMMON` lines can drift out of sync — a different order or type in one unit reinterprets the shared bytes with no error — whereas a single module declared once and `use`d everywhere has exactly one definition the compiler enforces across all users.

What's Next

You have modernized PLATE by reasoning through the eight steps. Chapter 19 turns that reasoning into a reference: a translation dictionary that tabulates every legacy pattern — COMMON, EQUIVALENCE, DATA, GOTO, arithmetic IF, statement functions, assumed-size arrays, implicit typing — beside its modern replacement, so that next time you meet old code you can look up the move instead of deriving it. It is the chapter you will keep open on a second monitor for the rest of your career, and it finishes the PLATE side quest using the dictionary to complete the steps your Project Checkpoint began.