Fortran subroutines and functions

As we start writing more complicated software to solve more complicated problems we will need to stay organized. The first tool to keep our larger projects tidy is the ability to split off commonly used code into functions and subroutines. This achieves two goals simultaneously:

  • This chunks potentially complicated code into smaller, more digestible pieces. Code written this way is much easier to read and understand (and debug!).

  • This also means that code that needs to be re-used does not need to be re-written for each use. You can easily just call the function or subroutine from wherever it is needed.

    • The upshot of this is that if you revise how a piece of code works you don’t need to go looking for every instance of it. Instead updating just the function/subroutine definition will accomplish your goal.

Functions, external and internal

Let’s dive straight into some examples, and from there look at the differences between functions and subroutines. Here’s an example of how to define and use your own function:

 1! /codes/fcn1.f90
 2
 3program fcn1
 4  implicit none
 5  real (kind=8) :: y, z
 6  real (kind=8), external :: f
 7
 8  y = 2.d0
 9  z = f(y)
10  print *, "y = ", y
11  print *, "z = ", z
12end program fcn1
13
14real (kind=8) function f(x)
15  implicit none
16  real (kind=8), intent(in) :: x
17  f = x**2
18end function f

Download this code

It prints out

y =    2.00000000000000
z =    4.00000000000000

Comments:

  • Functions return only one value, and by default the name of the return value is the same as the function name, see line 17

    • The name of the return value can be different from the function name, which we will see momentarily

  • f is declared external in the main program to let the compiler know it is a function defined elsewhere rather than a variable

  • Before the function keyword on line 14 we need to declare the return type

    • This will change when we specify the return variable explicitly

  • We need the ever-present implicit none statement to cover the scope of the function

  • After implicit none we need to declare the types for all of the input variables

    • You could move implicit none after the argument declarations, this is mostly a matter of style as long as you declare all of your arguments

  • The intent(in) statement tells the compiler that x is a variable passed into the function that will not be modified in the function

  • Here the program and function are in the same file. Later, we will see how to break things up so that functions and subroutines can be defined in another file

Alternatively, we can place the function inside our main program like this:

 1! /codes/fcn1_internal.f90
 2
 3program fcn1
 4  implicit none
 5  integer, parameter :: fp = selected_real_kind(15)
 6  real (fp) :: y,z
 7
 8  y = 2.0_fp
 9  z = f(y)
10  print *, "y = ", y
11  print *, "z = ", z
12
13contains
14
15  real (fp) function f(x)
16    implicit none
17    real (fp), intent(in) :: x
18    f = x**2
19  end function f
20
21end program fcn1

Download this code

Comments:

  • Our function is no longer declared external

  • Our function is now defined after it gets used

  • For simple functions like this one, this pattern is preferred over the previous version. Better yet, you should consider putting function definitions away in their own module (more on this below).

  • Why is it that we could use fp to specify our precision in this case, but needed kind=8 in the previous one?

As a final iteration let’s specify the return variable explcitly. This can be done as follows:

 1! /codes/fcn1_internal_retval.f90
 2
 3program fcn1
 4  implicit none
 5  integer, parameter :: fp = selected_real_kind(15)
 6  real (fp) :: y,z
 7
 8  y = 2.0_fp
 9  z = square_number(y)
10  print *, "y = ", y
11  print *, "z = ", z
12
13contains
14
15  ! Function that returns the square of its input
16  function square_number(x) result(sqx)
17    implicit none
18    real (fp), intent(in) :: x ! Number to square
19    real (fp) :: sqx           ! Return variable
20    sqx = x**2
21  end function square_number
22
23end program fcn1

Download this code

Comments:

  • The return type has moved into the function body next to the input type declarations

  • We can now give the function a nice clear name, but use a simpler variable inside for readability

  • We have also documented this function by providing comments desribing its purpose and all of its arguments

    • This is a very good habit to get into

    • You can even write structured comments that can aid documentation and editor features

Overall, this final variant is the one you should use for all function definitions. It is substantially easier to read and understand.

Modifying arguments

The input argument(s) to a function might also be modified, though this is not very common and is arguably bad practice (you should use a subroutine). Regardless, here is an example:

 1! /codes/fcn2.f90
 2
 3program fcn2
 4  implicit none
 5  real (kind=8) :: y, z
 6  real (kind=8), external :: f
 7
 8  y = 2.
 9  print *, "Before calling f: y = ", y
10  z = f(y)
11  print *, "After calling f:  y = ", y
12  print *, "z = ", z
13end program fcn2
14
15real (kind=8) function f(x)
16  implicit none
17  real (kind=8), intent(inout) :: x
18  f = x**2
19  x = 5.
20end function f

Download this code

This produces

Before calling f: y =    2.00000000000000
After calling f:  y =    5.00000000000000
z =    4.00000000000000

The use of intent

You are not required to specify the intent of each argument, but there are several good reasons for doing so:

  • It makes your code more readable to other people. Immediately seeing the intent of each argument makes it much easier to reason about someone else’s code.

  • It helps catch bugs. If you specify intent(in), and then the variable is changed in the function or subroutine, the compiler will give an error.

  • It can help the compiler produce machine code (e.g., *.o files) that runs faster. For example, if the optimizer knows that some variables will never change during execution, this fact can be used.

Subroutines

Functions are expected to produce a single output variable and examples like the one just given where an argument is modified are considered bad programming style.

Subroutines are more flexible since they can have any number of inputs and outputs, or might not have any output variables whatsoever. For example, a subroutine might take an array as an argument and store the array to some file on disk without returning anything to the calling program.

Here is a version of the same program as fcn1 above, which uses a subroutine instead:

 1! /codes/sub1.f90
 2
 3program sub1
 4  implicit none
 5  real (kind=8) :: y, z
 6
 7  y = 2.
 8  call fsub(y,z)
 9  print *, "z = ",z
10end program sub1
11
12subroutine fsub(x,f)
13  implicit none
14  real (kind=8), intent(in) :: x
15  real (kind=8), intent(out) :: f
16  f = x**2
17end subroutine fsub

Download this code

Comments:

  • Now we tell the compiler that x is an input variable and y is an output variable for the subroutine.

  • The intent declarations are optional but strongly recommended (see the above discussion)

  • Notice that fsub is not declared external, actually it isn’t declared at all prior to its use

  • What happens if we use contains to move the declaration inside our main program? What compiler flags could we use to test this change?

Here is a version that works on an array x instead of a single value:

 1! /codes/sub2.f90
 2
 3program sub2
 4  implicit none
 5  real (kind=8), dimension(3) :: y, z
 6  integer :: n
 7
 8  y = (/2., 3., 4./)
 9  n = size(y)
10  call fsub(n,y,z)
11  print *, "z = ",z
12end program sub2
13
14subroutine fsub(n,x,f)
15  ! compute f(x) = x**2 for all elements of the array x 
16  ! of length n.
17  implicit none
18  integer, intent(in) :: n
19  real (kind=8), dimension(n), intent(in) :: x
20  real (kind=8), dimension(n), intent(out) :: f
21  f = x**2
22end subroutine fsub

Download this code

This produces

4.00000000000000        9.00000000000000        16.0000000000000

Comments:

  • The length of the array is passed into the subroutine as an argument. Note that this length is used to specify the dimension of the remaining arguments.

    • In Fortran 77 this was the only way to pass arrays as arguments, so many codes continue to use this pattern to match any legacy calls.

  • In Fortran 90 (and beyond) we can write subroutines and functions with array type arguments without passing the array size. This is preferred in modern codes. (See the example below)

    • Doing this requires that the function or subroutine be defined in a module or placed inside the program using contains.

    • This has to do with implicit and explicit interfaces.

Subroutine in a module

Here is a version that avoids passing the length of the array into the subroutine. In order for this to work some additional interface information must be specified. This can be done by including the subroutine in a module.

 1! /codes/sub3.f90
 2
 3module sub3module
 4
 5contains 
 6
 7  subroutine fsub(x,f)
 8    ! compute f(x) = x**2 for all elements of the array x. 
 9    implicit none
10    real (kind=8), intent(in) :: x(:)
11    real (kind=8), intent(out) :: f(size(x))
12    f = x**2
13  end subroutine fsub
14
15end module sub3module
16
17!----------------------------------------------
18
19program sub3
20  use sub3module
21  implicit none
22  real (kind=8) :: y(3), z(3)
23
24  y = (/2., 3., 4./)
25  call fsub(y,z)
26  print *, "z = ",z
27end program sub3

Download this code

Comments:

  • See the section Fortran modules for more information about modules. Note, in particular, that the module must be defined first and then used in the program via the use statement.

    • Typically modules will be placed in their own source files, rather than in the same file as the main program.

  • The size intrinsic function lets you determine the array sizes internal to a subroutine, and makes it easy to write dimension independent code (but not rank independent).

  • Note that the use of the dimension keyword has been suppressed here compared to the previous example. Both syntaxes are valid, though you should choose one and stick to it for consistency.

  • The declaration of x in the subroutine uses x(:) indicates that it is an array with a single index (having rank 1), without specifying how long the array is. If x was a rank 2 array, indexed by x(i,j), then the declaration statement would use (:,:).

  • The declaration of f in the subroutine uses (size(x)) to indicate that it is an array with one index and the length of the array is the same as that of x. In fact it would be sufficient to tell the compiler that it is an array of rank 1.

    real (kind=8), dimension(:), intent(out) :: f
! Or
real (kind=8), intent(out) :: f(:)

    but indicating that it has the same size as x is useful for someone trying to understand the code, and may enable certain compiler optimizations.