Chapter 6: Procedures — Subroutines, Functions, and Organizing Your Code

"The purpose of abstraction is not to be vague, but to create a new semantic level in which one can be absolutely precise." — Edsger W. Dijkstra

Overview

Up to now your programs have been single, straight-line stories: declare some variables, loop over an array, print a result, stop. That is fine for a page of code. It is a disaster for the forty-thousand-line simulation you are working toward, because a program written as one long block cannot be tested in pieces, cannot be reused, and cannot be reasoned about by anyone — including you, six months from now. The cure, and the single most important organizing tool in the language, is the procedure: a named, self-contained unit of computation that you can call by name, hand some inputs, and get some outputs back.

Every serious Fortran code is a few thousand procedures wearing a trench coat. The weather model computes a radiation term by calling a radiation procedure; the linear-algebra library solves your system by calling dgesv; your own solver, by the end of this book, will march time forward by calling step. Procedures are how a computation is decomposed into pieces small enough to be correct, and then composed back into something large enough to be useful. This chapter is where your heat program stops being a script and starts being software.

Fortran's procedures also carry a feature that will quietly protect you for the rest of your career and that most languages simply do not have: intent, a declaration on every argument saying whether the procedure may read it, write it, or both — checked by the compiler, before your code ever runs. We will make a great deal of this, because it is a small habit with an enormous payoff in correctness, and because the same information that keeps you safe also helps the compiler make your code fast.

In this chapter, you will learn to:

• Decide, for any task, whether it wants a subroutine or a function, and write either one cleanly.
• Put intent(in), intent(out), or intent(inout) on every dummy argument — our universal house rule — and see the compiler reject a whole class of bugs because of it.
• Give procedures optional and keyword arguments, so one routine can serve simple and elaborate callers.
• Mark procedures pure and elemental when they qualify, and understand why that hands the optimizer more freedom (a theme we pay off in Chapter 27).
• Nest helper procedures inside a program with contains, and write a procedure that calls itself.
• Pass arrays the modern way — as assumed-shape arguments — and know why the old alternatives are traps.

Learning Paths

How to read this chapter by track. - 🔬 Scientist — §6.1, §6.2, and §6.6 are the load-bearing sections for your code; intent (§6.2) and assumed-shape arrays (§6.6) will prevent more bugs than any debugger. Read §6.4 for the payoff. - 📖 Standard — read straight through; procedures are the backbone of every later chapter. - 🔧 Legacy — §6.6's contrast between assumed-shape and the legacy assumed-size array, and the note on explicit interfaces in §6.5, are exactly the traps you will meet in old code (previewing Part IV). - ⚡ HPC — §6.4 (pure/elemental) is your section: these attributes are what let the compiler vectorize and parallelize calls. Skim §6.1; you know what a function is.

6.1 Subroutines and Functions: Two Tools, One Idea

A procedure is a named block of code you can invoke from elsewhere. Fortran gives you two kinds, and the difference between them is not arbitrary — it maps onto a real distinction in what you are asking the computer to do.

A function computes and returns a single value, and you use it inside an expression, exactly the way you use sqrt or sin:

area = pi * radius**2
energy = kinetic(mass, velocity) + potential(mass, height)

A subroutine performs an action and communicates through its arguments; you invoke it with the call statement, and it may return zero, one, or many results:

call solve(matrix, rhs, solution)      ! fills in `solution`
call write_field(temperature, "out.dat")

The rule of thumb is almost grammatical. A function is a noun — it is a value (the mean of this array, the area of that circle) and belongs in an expression. A subroutine is a verb — it does something (solve this system, write this file, update that field) and stands on its own as a statement. When a task produces one clean value and no side effects, reach for a function. When it produces several results, modifies its inputs in place, or performs an action like I/O, reach for a subroutine.

💡 Intuition: ask yourself, "could I write this on the right-hand side of an =?" x = mean(v) reads naturally, so mean is a function. x = solve(a, b, sol) reads like nonsense — solving fills in sol, it is not a value you assign — so solve is a subroutine.

A vocabulary note we will lean on for the rest of the book: the names a procedure declares for its inputs are its dummy arguments, and the values a caller actually supplies are the actual arguments. At a call, Fortran associates them by position — the first actual with the first dummy, the second with the second, and so on — unless you override that with the keyword syntax of §6.3. A dummy argument is not a free scratch variable; it is the compiler's stand-in for the caller's data, which is exactly why declaring what you may do to it — its intent (§6.2) — carries so much weight.

Here is the shape of both, in one small program that computes descriptive statistics. Notice the contains statement: everything after it, up to end program, is a procedure that belongs to the program. We will return to contains properly in §6.5; for now, read it as "and here are the helpers this program uses."

program stats_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: sample(5) = [2.0_dp, 4.0_dp, 4.0_dp, 4.0_dp, 6.0_dp]
  real(dp) :: mu, sigma

  ! A FUNCTION returns one value; use it directly in an expression.
  print '(a, f6.3)', 'mean (function)      = ', mean(sample)

  ! A SUBROUTINE returns several values through its arguments; call it.
  call describe(sample, mu, sigma)
  print '(a, f6.3)', 'mean (subroutine)    = ', mu
  print '(a, f6.3)', 'std dev (subroutine) = ', sigma

contains

  function mean(x) result(m)
    real(dp), intent(in) :: x(:)     ! read-only (§6.2); assumed-shape (§6.6)
    real(dp) :: m
    m = sum(x) / real(size(x), dp)
  end function mean

  subroutine describe(x, avg, sd)
    real(dp), intent(in)  :: x(:)
    real(dp), intent(out) :: avg, sd
    avg = mean(x)
    sd  = sqrt(sum((x - avg)**2) / real(size(x), dp))
  end subroutine describe

end program stats_demo

$ gfortran -std=f2018 -Wall -O2 example-01-subroutine-vs-function.f90 -o stats && ./stats
mean (function)      =  4.000
mean (subroutine)    =  4.000
std dev (subroutine) =  1.265

Three details of style, each of which we hold to for the rest of the book. First, the result clause: function mean(x) result(m) names the variable that carries the return value. Fortran also lets you assign to the function's own name instead, but the result form is clearer, and it becomes mandatory for the recursive functions of §6.5, so we use it everywhere. Second, mean returns one number, so it is a function; describe returns two (a mean and a standard deviation), so it is a subroutine — Fortran functions return a single result, and this is the everyday reason you choose a subroutine. Third, look at intent(in) and intent(out) on every argument. That is not decoration, and it is the subject of the next section.

One principle about functions is worth taking seriously from the start: a function should compute its value and nothing else. Fortran will technically permit a function to modify its arguments or print to the screen, but a function you cannot call twice without changing the program's behavior betrays every reader who reasonably expects f(x) to be just a value. Keep functions clean and push genuine actions into subroutines, where a reader expects side effects to live. (Section 6.4 gives you a way to promise a function is clean — the pure attribute — and have the compiler hold you to it.) You will also meet, in existing code, the older style that assigns to the function's own name in place of a result variable; that form is perfectly legal, but we prefer result for its clarity, and recursion (§6.5) requires it.

🐍 Python Comparison: In Python you would return several values as a tuple — return mean, std — and unpack them with mu, sigma = describe(x). Fortran functions return exactly one object, so multiple outputs travel through subroutine arguments instead. It feels heavier at first, but the argument list becomes an explicit, compiler-checked contract: each output is named, typed, and marked with its direction of data flow. Python trades that checking for brevity; in numerical code that runs for days, the checking is usually the better trade.

From History. Functions and subroutines were not in the very first FORTRAN — they arrived in FORTRAN II (1958), which introduced separately compiled subprograms. That single addition is what let programs grow past what one person could hold in their head, and it is why "modular" scientific codes have been possible for almost the entire history of the field.

6.2 intent: The Contract the Compiler Enforces

Here is a feature you will come to love, and one that C, C++, Python, Fortran-before-90, and most other languages you know cannot match. When you declare a dummy argument — the placeholder name a procedure uses for whatever the caller passes — you also declare its intent: whether the procedure will read it, write it, or both.

Definition (intent). An attribute on a dummy argument that states how the procedure uses it. intent(in) means the argument is read-only — the procedure may use its value but must not change it. intent(out) means write-only — the argument arrives undefined, and the procedure is expected to set it. intent(inout) means read and write — the procedure receives a meaningful value and may modify it in place. The compiler checks that the procedure's body honors whatever intent you declared.

Our house rule, stated in the style bible and followed in every example in this book, is simple: declare an intent on every dummy argument, always. It costs one word per argument and it buys you three things — better error-checking, clearer documentation, and better optimization — for that price.

The error-checking is the immediate payoff. If you declare an argument intent(in) and then, by accident, assign to it, the program does not compile. The mistake is caught at build time, at the exact line, with a clear message — not discovered three hours into a run as a mysteriously corrupted result.

🚪 Threshold Concept. intent turns a comment into a contract. In most languages, "this function does not modify its input" is a hope you write in a docstring and pray callers respect. In Fortran it is a declaration the compiler enforces on the function itself: an intent(in) argument physically cannot be modified, and the build fails if you try. Once you internalize that the argument list is a checked specification — not just a list of names — you start designing procedures by first deciding the direction of every piece of data. That habit, more than any single syntax rule, is what makes large Fortran codes maintainable.

The three intents in one place — a swap routine is the classic intent(inout), because it both reads and writes both arguments:

program swap_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: a, b
  a = 1.0_dp
  b = 2.0_dp
  call swap(a, b)
  print '(a, 2f6.2)', 'after swap: ', a, b

contains

  subroutine swap(x, y)
    real(dp), intent(inout) :: x, y    ! read AND write both
    real(dp) :: tmp
    tmp = x
    x   = y
    y   = tmp
  end subroutine swap

end program swap_demo

$ gfortran -std=f2018 -Wall -O2 swap_demo.f90 -o swap && ./swap
after swap:   2.00  1.00

Now watch the compiler earn its keep. Suppose we had carelessly marked swap's arguments intent(in):

subroutine swap_broken(x, y)
  real(dp), intent(in) :: x, y     ! WRONG: we promise not to change x, y ...
  real(dp) :: tmp
  tmp = x
  x   = y                          ! ... then change them anyway
  y   = tmp
end subroutine swap_broken

🐛 Find the Bug. The routine above will not compile. gfortran reports something like Error: Dummy argument 'x' with INTENT(IN) in variable definition context (assignment) at the line x = y. This is exactly the point of intent: you told the compiler x was read-only, and it is holding you to it. The fix is to declare the intent that matches what the code actually does — intent(inout) — or, if the modification was itself the mistake, to leave the intent and remove the assignment. Either way, a bug that in C would compile cleanly and corrupt data at run time is caught here before the program exists.

The same protection reaches the caller. Because swap promises to write its arguments, you cannot pass it something unwritable: call swap(1.0_dp, 2.0_dp) — swapping two literal constants — is a compile error, because a constant is not a "definable" thing. The compiler has caught a nonsensical call for free.

intent scales up to arrays and derived types unchanged, and there it earns even more. An intent(in) array argument is a read-only view of the caller's data — the procedure may compute from it but cannot alter it — which is one of the guarantees that lets the compiler assume a routine's input and output do not secretly overlap. An intent(out) array is fully re-defined by the procedure, and an intent(inout) array is updated in place, exactly as field is in the solver's step you will write at the end of this chapter. One wrinkle is worth knowing before you meet it: for a derived type with default component values (Chapter 9), an intent(out) argument is reset to those defaults on entry and any allocatable components are deallocated — a convenience once you are expecting it, and a genuine surprise if you are not.

There is one sharp edge, and it is worth meeting now rather than in a debugger.

⚠️ Common Pitfall: an intent(out) argument arrives undefined. Whatever value the caller's variable held is gone the moment the procedure begins — Fortran is entitled to forget it. So if you write a routine meant to add to a running total but declare the accumulator intent(out), you silently discard the incoming total every call. The rule: use intent(out) only for a result you compute from scratch; if the procedure needs the incoming value, it is intent(inout). (For a plain number, reading an intent(out) argument before you set it is undefined behavior; the flag -fcheck=all, which we meet in Chapter 13, helps catch use-before-definition during development.)

The third benefit — optimization — is quieter but real. When the compiler knows an argument is intent(in), it knows the procedure will not write through it, which is one of the guarantees that lets Fortran generate famously fast numerical code. We foreshadowed this as the "no-aliasing advantage" in Chapter 1, and we will measure it in Chapter 27. For now, notice the pattern that will recur through the whole book: the same declaration that keeps you safe also makes you fast. intent is information, and in Fortran, information is speed.

🔄 Check Your Understanding 1. You are writing subroutine normalize(v, scale) that divides the vector v by scale in place and leaves scale untouched. What intent belongs on each argument? 2. Why does call swap(1.0_dp, 2.0_dp) fail to compile, while call swap(a, b) succeeds? 3. A colleague's routine takes an intent(out) array and only fills in its interior, expecting the edges to keep their previous values. Why is that a bug?

6.3 Optional and Keyword Arguments

Real procedures often want to serve two kinds of caller: the one who wants the simple default, and the one who wants to tune every knob. Fortran supports both with optional arguments — arguments the caller may supply or omit — and keyword arguments — the ability to name arguments at the call site instead of relying on their position.

Definition (optional argument). A dummy argument declared with the optional attribute, which the caller may leave out. Inside the procedure, the intrinsic present(arg) returns .true. if the caller supplied arg and .false. if they did not, so you can substitute a default. You must never use the value of an absent optional argument — always guard it with present first.

Definition (keyword argument). At a call, naming a dummy argument explicitly with name=value, e.g. call relax(x, old, goal, factor=0.5_dp). Keyword arguments may be given in any order and make it possible to supply a later optional argument while skipping an earlier one. Once you use a keyword in a call, every argument after it must also use a keyword.

A worked example ties both together. This relax routine nudges a value a fraction of the way toward a target — the essential move of every iterative solver, including the heat solver you are building. The fraction, factor, is optional and defaults to a full step of 1.0:

program optional_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: a, b, c

  call relax(a, 100.0_dp, 0.0_dp)                          ! factor omitted -> 1.0
  call relax(b, 100.0_dp, 0.0_dp, 0.5_dp)                  ! factor positional
  call relax(c, x_old=100.0_dp, goal=0.0_dp, factor=0.1_dp)  ! all by keyword

  print '(a, f7.2)', 'full step  (factor 1.0) : ', a
  print '(a, f7.2)', 'half step  (factor 0.5) : ', b
  print '(a, f7.2)', 'tenth step (factor 0.1) : ', c

contains

  subroutine relax(x_new, x_old, goal, factor)
    real(dp), intent(out)          :: x_new
    real(dp), intent(in)           :: x_old, goal
    real(dp), intent(in), optional :: factor
    real(dp) :: f

    f = 1.0_dp                      ! the default ...
    if (present(factor)) f = factor ! ... overridden only if supplied
    x_new = x_old + f * (goal - x_old)
  end subroutine relax

end program optional_demo

$ gfortran -std=f2018 -Wall -O2 example-02-intent-and-optional.f90 -o relax && ./relax
full step  (factor 1.0) :    0.00
half step  (factor 0.5) :   50.00
tenth step (factor 0.1) :   90.00

Read the three calls carefully, because they show the whole feature. The first omits factor entirely, so present(factor) is .false. and f keeps its default of 1.0 — a full step from 100 all the way to the goal of 0. The second passes 0.5 positionally, so we relax halfway, landing at 50. The third names every argument, which lets us pass them in any order and documents the call at the point of use — goal=0.0_dp is far clearer than a bare 0.0_dp third in a list. All three go through the same routine.

🐍 Python Comparison: this is Fortran's version of Python's def relax(x_old, goal, factor=1.0). The differences are instructive. Python attaches the default in the signature; Fortran has no default-value syntax, so you write the default in the body and gate it with present. Python lets any argument be passed by keyword; Fortran does too, and with the same rule that positional arguments must come before keyword ones. The one thing Fortran can do that Python cannot is distinguish "passed the default value" from "did not pass it at all" — present tells you which, which matters when the default is expensive to compute or when "absent" should mean something different from any value.

⚠️ Common Pitfall: never touch an absent optional's value. Writing x_new = x_old + factor * ... without the present guard reads factor even when the caller omitted it — undefined behavior, and a classic source of nondeterministic bugs. The discipline is unconditional: for every optional argument, either test present before using it, or pass it straight through to another procedure that will (an absent optional may be forwarded as an absent optional). If you find yourself wanting the value with no default in mind, the argument probably should not have been optional.

Optional arguments are not only for input defaults. A common and genuinely useful pattern is the optional output — a diagnostic the caller may or may not want. A linear solver might declare real(dp), intent(out), optional :: residual and compute it only when asked, guarding the assignment with if (present(residual)). The caller who cares passes a variable and reads the residual back; the caller who does not omits it and pays nothing for a number they were going to ignore. And when a wrapper procedure needs to relay an optional argument to an inner one, Fortran lets you pass it along even when it is absent — "absent" is forwarded as "absent" — so the wrapper need not test present itself, as long as the inner routine also declares that argument optional. This is how thin convenience wrappers stay thin.

6.4 pure and elemental: Promises That Make Code Fast

Some procedures are mathematically clean: hand them the same inputs and they always return the same outputs, and they touch nothing else in the program along the way — no file writes, no global variables, no surprises. Fortran lets you declare that cleanliness, and when you do, the compiler can optimize far more aggressively because it knows the call has no hidden consequences.

Definition (pure). A procedure marked pure promises it has no side effects: it does not perform I/O, does not modify any global or module variable, and modifies nothing outside its own outputs. Every dummy argument of a pure function must be intent(in). The compiler enforces all of this — try to print inside a pure procedure and the build fails. Because a pure call cannot affect anything but its result, the compiler may reorder such calls, hoist them out of loops, or evaluate them in parallel.

Definition (elemental). A procedure marked elemental is written for scalar arguments but may be called with array arguments of any shape, in which case it is applied independently to each element. An elemental procedure is automatically pure. This is exactly how the intrinsics behave — sqrt(x) works whether x is a scalar or a million-element array — and elemental lets you write your own functions that do the same.

The two ideas in one program. to_kelvin is elemental, so the same function converts a single temperature or a whole array of them, with no loop on your part. average is pure and takes an assumed-shape 2-D array (the subject of §6.6):

program pure_elemental_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: temps_c(4) = [-40.0_dp, 0.0_dp, 37.0_dp, 100.0_dp]
  real(dp) :: temps_k(4)
  real(dp) :: grid(2, 3)
  integer  :: i, j

  ! elemental: one definition, works on a scalar AND on a whole array
  print '(a, f8.2)',  'body temp in K    : ', to_kelvin(37.0_dp)
  temps_k = to_kelvin(temps_c)               ! applied element by element, no loop
  print '(a, 4f8.2)', 'Celsius -> Kelvin : ', temps_k

  do j = 1, 3
     do i = 1, 2
        grid(i, j) = real(i + j, dp)
     end do
  end do
  print '(a, f8.2)',  'grid average      : ', average(grid)

contains

  elemental function to_kelvin(celsius) result(kelvin)
    real(dp), intent(in) :: celsius        ! scalar only — that is the rule
    real(dp) :: kelvin
    kelvin = celsius + 273.15_dp
  end function to_kelvin

  pure function average(a) result(avg)
    real(dp), intent(in) :: a(:,:)         ! assumed-shape: any 2-D array
    real(dp) :: avg
    avg = sum(a) / real(size(a), dp)
  end function average

end program pure_elemental_demo

$ gfortran -std=f2018 -Wall -O2 example-03-pure-elemental-assumed-shape.f90 -o pe && ./pe
body temp in K    :   310.15
Celsius -> Kelvin :   233.15  273.15  310.15  373.15
grid average      :     3.50

The elemental function is the star here. Written once for a single real, it applies to the four-element temps_c array with no loop, no size, no indexing — temps_k = to_kelvin(temps_c) reads like mathematics and compiles to a loop the compiler is free to vectorize. That is the same whole-array thinking you met in Chapter 5, now extended to your own functions. The one constraint to remember: an elemental procedure's arguments must be scalars. You are describing what to do to one element; the language handles the fanning-out across an array.

⚡ Performance Note: why does the optimizer care? A pure function is a promise that calling it changes nothing but its result. That promise lets the compiler do things it otherwise could not: evaluate the call once instead of twice if it appears twice, move it out of a loop if its arguments do not change, and — crucially — run it on many array elements at the same time, because independent elements cannot interfere. This is precisely the property that lets a do concurrent loop (Chapter 4) or an OpenMP-parallelized loop (Chapter 33) run safely across cores. pure and elemental are not cosmetic labels; they are optimization licenses you grant the compiler, and Chapter 27 shows what it does with them.

Make that freedom concrete. Imagine a loop whose body reads y(i) = to_kelvin(t) * a(i) + to_kelvin(t) * b(i), with t unchanged across the loop. Because to_kelvin is pure, the compiler knows both calls must return the same value and that calling it disturbs nothing else, so it may legally compute to_kelvin(t) once, before the loop, and reuse the result — an optimization it could not risk if the function might, say, advance a random-number generator or write a log line on each call. That is the whole game: pure removes the compiler's doubt. Every guarantee you can make about a procedure is a degree of freedom you hand the optimizer, and in numerical code those degrees of freedom convert directly into speed.

A practical habit worth forming now: make your small mathematical helpers pure, and your scalar-to-scalar transformations elemental, whenever they honestly qualify. It costs nothing, documents the routine's cleanliness, and keeps the door open to parallelism you may add later. The one caution is honesty — the compiler will reject a pure procedure that does I/O or touches global state, and it is right to. If a routine genuinely needs a side effect, it is not pure, and forcing the label is a mistake, not an optimization.

6.5 Internal Procedures and Recursion

Where should a procedure live? You have already been using the answer: after a contains statement, inside the program itself.

Definition (internal procedure). A procedure defined inside another program unit — a program, or (from Chapter 8) a module — after its contains statement. An internal procedure can see the host's variables directly, a feature called host association, and — importantly — the host always knows the internal procedure's exact interface.

That last point is not a technicality; it is the reason we put procedures inside contains at all. When a procedure is internal, the compiler has its full explicit interface on hand at every call: it knows the number, type, and intent of every argument, and it checks each call against them. Pass the wrong number of arguments, or a real where an integer is expected, and the build fails. The features of this chapter that feel most modern — optional and keyword arguments (§6.3), and assumed-shape arrays (§6.6) — require an explicit interface to work at all, and contains provides one automatically.

🔗 Connection: the historical alternative is the external procedure: a standalone subprogram with no host, whose interface the caller cannot see. Mismatched calls to external procedures — a wrong argument type, a missing argument — go uncaught and produce garbage at run time, and this is one of the great hazards of old FORTRAN 77 code you will meet in Part IV. The modern fix is exactly the one this chapter teaches: put procedures where they get an explicit interface — inside contains, or, better for anything reused across files, inside a module, which Chapter 8 is entirely about. For now, contains inside your program is enough.

Host association is powerful and occasionally dangerous. An internal procedure can read and write the host's variables without them being passed as arguments — convenient for a quick helper, but a trap if you forget and modify a host variable by accident. The safe default is to pass what a procedure needs as arguments, with explicit intents, and lean on host association sparingly for genuinely shared context. Note also that internal procedures do not nest: a procedure inside contains cannot itself have a contains of its own.

Recursion. A procedure may call itself. The canonical illustration is the factorial, $n! = n \times (n-1)!$, with the base case $0! = 1! = 1$:

🧩 Try It Yourself: predict this program's output before you compile it. The factorial grows fast, so we return a 64-bit integer to leave room.

```fortran program factorial_demo use, intrinsic :: iso_fortran_env, only: int64 implicit none integer :: n do n = 0, 6 print '(a, i1, a, i4)', 'factorial(', n, ') = ', factorial(n) end do

contains

recursive function factorial(n) result(f) integer, intent(in) :: n integer(int64) :: f if (n <= 1) then f = 1_int64 ! base case else f = int(n, int64) * factorial(n - 1) ! recursive case end if end function factorial

end program factorial_demo ```

console $ gfortran -std=f2018 -Wall -O2 factorial_demo.f90 -o fact && ./fact factorial(0) = 1 factorial(1) = 1 factorial(2) = 2 factorial(3) = 6 factorial(4) = 24 factorial(5) = 120 factorial(6) = 720

Two points of syntax and one of judgment. Syntactically, note the recursive prefix, and note that a recursive function must use the result clause — without it, the name factorial inside the body would be ambiguous between "the return value" and "call myself again." (Fortran 2018 actually makes procedures recursive by default, so the recursive keyword is technically optional under our -std=f2018 baseline; we write it anyway, because it is clear, and because it compiles under every standard.)

The judgment is a performance one, and it is honest. Recursion is elegant, but for the tight numerical kernels this book is about, iteration is almost always the right choice. A loop has no per-call overhead and no risk of exhausting the call stack; a deeply recursive routine can overflow the stack and crash. Factorial is a fine teaching example and a poor production one — you would compute it with a do loop from Chapter 4. Recursion earns its place for genuinely tree-shaped problems (traversing a hierarchy, some divide-and-conquer algorithms), not for things a loop does cleanly. Reach for it when the problem is recursive, not merely because the language allows it.

A fair example of recursion pulling its weight is a binary search over a sorted array: each call inspects the middle element and then recurses into just the left or right half, so the problem genuinely shrinks into a smaller copy of itself, and the recursion depth is only about $\log_2 n$ — a couple of dozen levels even for a billion elements, nowhere near a stack overflow. Contrast summing an array, where the "smaller copy" is merely one element shorter and the depth would equal the array's length: a loop is plainly right. The test is the shape of the problem, not the size of the data.

🔄 Check Your Understanding 1. Name two things a procedure gets by being internal (inside contains) that an external procedure does not. 2. Why must a recursive function use a result clause? 3. Give one reason to prefer an iterative loop over recursion in a numerical kernel.

6.6 Array Arguments: Assumed-Shape and Its Rivals

Passing arrays to procedures is where Fortran's design pays off most visibly — and where the language has accumulated three different mechanisms across sixty years, only one of which you should reach for by default. Getting this right is the difference between array code that is safe and self-describing and array code that silently corrupts memory.

Definition (assumed-shape array). A dummy array argument declared with only colons for its dimensions — real(dp), intent(in) :: a(:,:) — whose shape is taken automatically from the actual argument at each call. Inside the procedure, size, shape, lbound, and ubound report the real extents, and whole-array operations and array sections work exactly as they do on a locally declared array. Assumed-shape arguments require an explicit interface (so the procedure must be internal or in a module), and the dummy's lower bounds default to 1 regardless of the caller's bounds.

The modern default is assumed-shape, and you have already been using it: mean(x(:)) in §6.1, average(a(:,:)) in §6.4. The shape rides along with the array, so you never pass dimensions by hand, and every array intrinsic works inside the procedure. Here is a matrix trace — the sum of the diagonal — written the modern way, discovering the matrix's size for itself:

program trace_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  ! reshape fills COLUMN-MAJOR (Chapter 5), so this is diag(1, 2, 3):
  real(dp) :: m(3,3) = reshape([1.0_dp, 0.0_dp, 0.0_dp, &
                                0.0_dp, 2.0_dp, 0.0_dp, &
                                0.0_dp, 0.0_dp, 3.0_dp], [3, 3])
  print '(a, f6.2)', 'trace = ', trace(m)

contains

  pure function trace(a) result(t)
    real(dp), intent(in) :: a(:,:)               ! assumed-shape
    real(dp) :: t
    integer  :: i
    t = 0.0_dp
    do i = 1, min(size(a, 1), size(a, 2))
       t = t + a(i, i)
    end do
  end function trace

end program trace_demo

$ gfortran -std=f2018 -Wall -O2 trace_demo.f90 -o trace && ./trace
trace =   6.00

One assumed-shape subtlety catches everyone exactly once. The dummy's lower bounds default to 1, no matter what bounds the actual argument had. If the caller declares real(dp) :: g(0:9) and passes it to a routine whose dummy is x(:), then inside that routine x runs from 1 to 10, not 0 to 9 — the shape (ten elements) is preserved, but the numbering restarts at 1. This is almost always what you want, and it is why solver code indexes from 1 by habit. When the mathematics genuinely wants a different origin — a stencil naturally indexed from 0, say — you can ask for it explicitly with x(0:), which keeps the extent from the caller but starts the numbering at 0. Reach for that only when the problem asks; the default of 1 is the path of least surprise.

Now meet the two older mechanisms, so you can recognize them and know why you are not using them.

Explicit-shape arrays declare their extents outright, which means the caller must pass the sizes as separate arguments:

pure function trace_explicit(a, n) result(t)
  integer,  intent(in) :: n
  real(dp), intent(in) :: a(n, n)     ! extents passed in by hand
  real(dp) :: t
  integer  :: i
  t = 0.0_dp
  do i = 1, n
     t = t + a(i, i)
  end do
end function trace_explicit

This works, and it has a place — an explicit-shape array is guaranteed contiguous, which occasionally matters for the last drop of performance (a topic for Chapter 29). But you must thread the dimensions through every call yourself, and every hand-carried n is a chance to pass the wrong one and read past the end of the array. The assumed-shape version simply cannot make that mistake, because it never asks you for n.

The third mechanism is a genuine relic, and you should know it only well enough to replace it on sight.

🔧 Modern vs Legacy: the FORTRAN 77 assumed-size array writes its last dimension as *, meaning "some length I will not tell you." The procedure literally cannot know the array's true extent — size does not work, whole-array operations do not work, and bounds checking is defeated.

fortran ! Legacy (FORTRAN 77): assumed-SIZE. No kind, no intent, size unknown. subroutine scale_old(a, n, factor) integer :: n, i real :: a(*), factor ! the '*' hides the real length do i = 1, n a(i) = a(i) * factor end do end subroutine scale_old

fortran ! Modern: assumed-SHAPE. The shape travels with the array; no separate n. subroutine scale_new(a, factor) real(dp), intent(inout) :: a(:) real(dp), intent(in) :: factor a = a * factor ! whole-array op — the shape is known end subroutine scale_new

The modern routine is shorter, needs no n, uses a whole-array operation, gets bounds checking under -fcheck=all, and carries explicit intents. The legacy one asks you to pass n correctly on faith and punishes you with silent corruption if you slip. When you modernize old code in Part IV, turning a(*) into a(:) is step five of the eight-step recipe.

So why is assumed-shape "almost always right"? Because it is safe (the compiler and -fcheck=all know the real bounds), self-describing (the size travels with the data, so there is nothing to pass wrong), expressive (every array intrinsic and whole-array operation works inside the procedure), and general (one routine handles any size of plate). The price is small and worth naming honestly.

⚡ Performance Note: an assumed-shape argument is passed as a small array descriptor — a handful of words recording the base address, extents, and strides — rather than a bare pointer. That indirection is normally negligible. It matters only in two situations: extremely hot inner loops where the tiny per-access overhead adds up, and cases where the actual argument is a non-contiguous section (say, a(:, 2:100:2)), which can defeat some vectorization. The escape hatches — the contiguous attribute (Chapter 11) and, rarely, explicit-shape arrays — exist for exactly these cases, and you deploy them only when a profiler (Chapter 28) tells you to, never preemptively. Write assumed-shape first; optimize the descriptor away only if measurement demands it.

🔄 Check Your Understanding 1. Inside a procedure with real(dp), intent(in) :: a(:,:), how do you find the number of rows? 2. What can an assumed-shape array do that a legacy assumed-size array (a(*)) cannot? 3. Why does an assumed-shape argument require the procedure to be internal or in a module?

Project Checkpoint

Your heat program currently updates the temperature field with a block of code sitting in the main program — the whole-array Laplacian you wrote in Chapter 5. That was fine while the program was a script. It is time to make it software: we lift the update out of main and into a named procedure, step, with an intent on every argument and an assumed-shape array — the interface every later chapter will build on.

The canonical signature, fixed for the rest of the book, is step(field, alpha, dt): field is the 2-D temperature array (assumed-shape, so the routine works for any plate size), alpha is the diffusivity, and dt is the time step. field is intent(inout) — we read the current temperatures and overwrite them; alpha and dt are intent(in).

subroutine step(field, alpha, dt)
  real(dp), intent(inout) :: field(:,:)     ! assumed-shape: any plate size
  real(dp), intent(in)    :: alpha, dt      ! read-only parameters
  real(dp), allocatable   :: old(:,:)
  real(dp) :: lap
  integer  :: i, j, nx, ny

  nx = size(field, 1)
  ny = size(field, 2)
  allocate(old(nx, ny))
  old = field                               ! update every cell from one snapshot

  ! PLACEHOLDER physics — the real five-point stencil, the CFL limit on dt, and
  ! the boundary conditions all arrive in Chapter 24. Here each interior cell
  ! just moves toward its four neighbours, so we can lock in the INTERFACE now.
  do j = 2, ny - 1
     do i = 2, nx - 1
        lap = old(i-1, j) + old(i+1, j) + old(i, j-1) + old(i, j+1) &
              - 4.0_dp * old(i, j)
        field(i, j) = old(i, j) + alpha * dt * lap
     end do
  end do
end subroutine step

Three deliberate choices, each a lesson from this chapter. It is a subroutine, not a function, because it updates the field in place and returns nothing (§6.1). Every argument carries an intent (§6.2), so the compiler guarantees step cannot accidentally scribble on alpha or dt, and a reader knows at a glance that only field comes back changed. And field is assumed-shape (§6.6), so the identical routine serves a $4\times4$ test plate and a $4000\times4000$ production run — you never pass the grid dimensions by hand. The full driver in code/project-checkpoint.f90 runs it on a tiny plate held hot along the top edge; after one step, heat has diffused exactly one row inward:

after one step:
  100.00  100.00  100.00  100.00
    0.00   10.00   10.00    0.00
    0.00    0.00    0.00    0.00
    0.00    0.00    0.00    0.00

Do not read too much physics into the placeholder — it is a stand-in whose only job is to make the interface concrete and compilable. The genuine stencil, the stability condition that bounds dt, and the boundary handling are the whole point of Chapter 24, and when you reach it, the signature of step will not change — only its body. That is the payoff of settling the interface now: every future chapter (modules in Ch. 8, a field_t type in Ch. 9, OpenMP in Ch. 33) drops into a step whose shape the rest of the code already trusts.

Summary

This chapter turned your program into a set of cooperating procedures and gave every argument a compiler-checked direction of flow.

The two things to memorize. First, an intent on every argument, always — it is the cheapest bug prevention in the language and the habit that most separates robust Fortran from fragile Fortran. Second, arrays go in as a(:,:) — assumed-shape is safe, self-describing, and works with every array intrinsic; the alternatives are for special cases you will recognize when you meet them.

Spaced Review

Three questions reaching back to Chapter 4 (control flow) and Chapter 5 (arrays), the ground this chapter's procedures stand on.

1. (Ch. 5) Inside step we wrote allocate(old(nx, ny)) and then old = field. What kind of array is old, what does allocate do, and when is its memory released?
   

   Answer

   `old` is an allocatable array. `allocate` reserves storage for it at run time with the requested shape (here matching `field`). Because it is a local allocatable, it is automatically deallocated when the procedure returns — no explicit `deallocate` is needed.

2. (Ch. 5) step uses size(field, 1) and size(field, 2). In Fortran's column-major layout, which of these counts rows, and why does looping with j (columns) outermost and i (rows) innermost match the memory order?
   

   Answer

   `size(field, 1)` is the first dimension — the rows. Fortran stores arrays column-major, so the first index varies fastest through memory; walking the inner loop over `i` (down a column) touches consecutive memory addresses, which is cache-friendly. The full performance story is [Chapter 27](../../part-07-performance/chapter-27-why-fortran-is-fast/index.md).

3. (Ch. 4) The relax routine chose its default with if (present(factor)) f = factor. Rewrite the decision using an if ... else ... end if block, and say when the one-line form is preferable.
   

   Answer

   if (present(factor)) then; f = factor; else; f = 1.0_dp; end if. The one-line logical `if` is preferable when there is a single simple action and no `else` branch (here we had already set the default), keeping the common case terse.

4. (Ch. 4 & 5) Our placeholder loops do j = 2, ny-1 / do i = 2, nx-1 deliberately skip the first and last row and column. In heat-equation terms, what do those skipped cells represent, and which Chapter 4 construct would you use to reset them to fixed values each step?
   

   Answer

   The skipped edge cells are the boundary — the plate's fixed- temperature edges (Dirichlet conditions, formalized in Ch. 24). You would reset them with array-section assignments (Ch. 5), e.g. `field(1,:) = 100.0_dp`, or a masked `where` (Ch. 4) if the condition were value-dependent.

What's Next

You can now decompose a computation into procedures and hand data across their boundaries safely. But your heat program is still one file with all its procedures crammed into a single contains, and its parameters are still hard-coded. Chapter 7 fixes the second problem: it teaches Fortran's I/O — formatted output, list-directed reading, files, and the wonderfully practical namelist, which will let your solver read alpha, dt, and the grid size from a configuration file instead of a recompile. After that, Chapter 8 solves the first problem by moving these procedures out of the program and into modules — where, not coincidentally, they get the same explicit interface that contains gave them here, and the solver finally splits into the clean, reusable pieces a real scientific code is made of.