Basics of Fortran

Note This chapter of the lecture note has been partially extracted and modified from Prof. LeVeque’s (Univ. of Washington) online note on Fortran. Some parts of the chapter also has been adopted and modified from Prof. Garaud’s (AMS, UCSC) short notes on Fortran:

History

FORTRAN stands for FORmula TRANslator and was the first major high level language to catch on. The first compiler was written in 1954-57. Before this, programmers generally had to write programs in a low-level programming language called assembly.

Many version followed: Fortran II, III, IV. Fortran 66 followed a set of standards formulated in 1966.

See also

for brief histories.

Fortran 77

The standards established in 1977 lead to Fortran 77, or f77, and many codes are still in use that follow this standard.

Fortran 77 does not have all the features of newer versions and many things are done quite differently.

One feature of f77 is that lines of code have a very rigid structure. This was required in early versions of Fortran due to the fact that computer programs were written on punched cards. All statements must start in column 7 or beyond and no statement may extend beyond column 72. The first 6 columns are used for things like labels (numbers associated with particular statements). In f77 any line that starts with a ‘c’ in column 1 is a comment.

We will not use f77 in this class but if you need to work with Fortran in the future you may need to learn more about it because of all the legacy codes that still use it (see also f77-wikibooks).

Fortran 90/95

Dramatically new standards were introduced with Fortran 90, and these were improved in mostly minor ways in Fortran 95. There are newer Fortran 2003 and 2008 standards but few compilers implement these fully yet. See Wikipedia page on Fortran standards for more information.

For this class we will use the Fortran 90/95 standards, which we will refer to as Fortran 90 for brevity.

See also online documents

for more online tutorials on Fortran 90 and Fortran 2003.

Why Fortran?

Frequently, people ask what is the advantage of using Fortran as opposed to using other modern scientific languages, such as C/C++. One good explanation can be found here: FAQ.

Compilers

Unlike Python code (which is an interpreted language; we will study this later), a Fortran program (which is a compiled language) must pass through several stages before being executed (i.e., compilation). There are several different compilers that can turn Fortran code into an executable (or binary), as described more below.

In this class we will use gfortran, which is an open source compiler, part of the GNU Project. See http://gcc.gnu.org/fortran/ for more about gfortran.

There is an older compiler in this suite called g77 which compiles Fortran 77 code, but gfortran can also be used for Fortran 77 code and has replaced g77.

There are several commercial compilers which are better in some ways, in particular they sometimes do better optimization and produce faster running executables. They also may have better debugging tools built in. Some popular ones are the Intel and Portland Group compilers.

See a list of Fortran compilers-wiki.

File extensions

For the gfortran compiler the typical convention is to use the lower-case file extensions. In this case fixed format code should use .f, while free format code should have the extension .f90 or .f95. We will use .f90 in this course.

Another convention is to use the upper-case extensions such as .F for fixed format code, and .F90 or .F95 for free format code.

The difference between the lower-case and the upper-case extensions is that the latter has an extended feature that can be preprocessed in a C-like style. This can be particularly very useful if one develops a code project where Fortran and C/C++ routines are both used and Fortran routines are to be preprocessed. See a FLASH code example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
!!****if* source/physics/Hydro/HydroMain/unsplit/hy_uhd_avgState
!!
!! NAME
!!
!!  hy_uhd_avgState
!!
!! SYNOPSIS
!!
!!  hy_uhd_avgState( integer(IN) :: sweepDir,
!!                   real(IN)  :: VL(HY_VARINUM3),
!!                   real(IN)  :: VR(HY_VARINUM3),
!!                   real(OUT) :: Vavg(HY_VARINUM2) )
!!
!! DESCRIPTION
!!
!!  This routine computes proper average state values at each interface
!!  using a simple arithmatic average.
!!  The calculated average state values are used in the Roe and Lax-Friedrichs
!!  Riemann solvers.
!!
!! ARGUMENTS
!!
!!  sweepDir - sweep direction
!!  VL    -  a vector array for the left state  
!!            (DENS,VELX,VELY,VELZ,MAGX,MAGY,MAGZ,PRES + GAMC,GAME,EINT)
!!  VR    -  a vector array for the right state 
!!            (DENS,VELX,VELY,VELZ,MAGX,MAGY,MAGZ,PRES + GAMC,GAME,EINT)
!!  Vavg  -  a vector array for the computed average state
!!            (DENS,VELX,VELY,VELZ,MAGX,MAGY,MAGZ,PRES + GAMC,GAME)
!!
!!*** 
#include "Flash.h"
#include "UHD.h"

Subroutine hy_uhd_avgState(sweepDir,VL,VR,Vavg)

#if defined(FLASH_USM_MHD) || defined(FLASH_UGLM_MHD)
  use Hydro_data,           ONLY : hy_forceHydroLimit
#endif
  use hy_uhd_slopeLimiters, ONLY : signum

  implicit none

#include "Flash.h"
#include "UHD.h"
  !! Arguments type declaration -----------------
  integer, intent(IN) :: sweepDir
  real, dimension(HY_VARINUM3), intent(IN)  :: VL,VR
  real, dimension(HY_VARINUM2), intent(OUT) :: Vavg
  !! --------------------------------------------
  real :: sig

  !! Arithmetic averages
  Vavg(HY_DENS:HY_GAME) = .5*(VL(HY_DENS:HY_GAME)+VR(HY_DENS:HY_GAME))

#if defined(FLASH_USM_MHD) || defined(FLASH_UGLM_MHD)
  if (hy_forceHydroLimit) Vavg(HY_MAGX:HY_MAGZ) = 0.
#endif

  !! Use upwinding for game and gamc that are averaged along the
  !! entropy wave.
  sig = signum(Vavg(sweepDir+1))
  Vavg(HY_GAMC:HY_GAME) = 0.5*( (1.+sig)*VL(HY_GAMC:HY_GAME) &
                               +(1.-sig)*VR(HY_GAMC:HY_GAME) )

End Subroutine hy_uhd_avgState

Download this code

You see in the above example the lines with macros starting with # in the first column. They are the C-like preprocessor macros that are found in lines:

  • 32, 33, 37-39, 44, 45, and 56-58

Compiling, linking, and running a Fortran code

Suppose for example we have a Fortran file named demo1.f90. We can not run this directly the way we would do in MATLAB or Python where one does not require to compile the script for execution. Instead, it must be converted into object code, a version of the code that is in a machine language specific to the type of computer. This is done by the compiler.

Then a linker must be used to convert the object code into an executable (or binary) that can actually be executed.

This is broken into two steps because often large programs are split into many different *.f90 files. Each one can be compiled into a separate object file, which by default has the same name but with a .o extension (for example, from demo1.f90 the compiler would produce demo1.o). One may also want to call on library routines that have already been compiled and reside in some library. The linker combines all of these into a single executable.

For more details on the process, see for example:

For the simplest case of a self-contained program in one file, we can combine both stages in a single gfortran command, e.g.

$ gfortran demo1.f90

By default this will produce an executable named a.out for obscure historical reasons (it stands for assembler output, see wikipedia).

To run the code you would then type

$ ./a.out

Note we type ./a.out to indicate that we are executing a.out from the current directory (see Basic Unix/Linux Commands). There is an environment variable PATH that contains your search path, the set of directories that are searched whenever you type a command name at the Unix prompt. Often this is set so that the current directory is the first place searched, in which case you could just type a.out instead of ./a.out. However, it is generally considered bad practice to include the current directory in your search path because bad things can happen if you accidentally execute a file.

If you don’t like the name a.out you can specify an output name using the -o flag with the gfortran command. For example, if you like the Windows convention of using the extension .exe for executable files

$ gfortran demo1.f90 -o demo1.exe
$ ./demo1.exe

will also run the code.

Note that if you try one of the above commands, there will be no file demo1.o created. By default gfortran removes this file once the executable (or binary, e.g., demo1.exe, or a.out) is created.

Later, we will see that it is often useful to split up the compile and link steps, particularly if there are several files that need to be compiled and linked. We can do this using the -c flag to compile without linking

$ gfortran -c demo1.f90              # produces demo1.o
$ gfortran demo1.o -o demo1.exe      # produces demo1.exe

There are many other compiler flags that can be used, see linux man page for gfortran for a list.

Sample codes

The first example simply assigns some numbers to variables and then prints them out. The comments below the code explain some features.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
! /codes/demo1.f90

program demo1

  ! Fortran 90 program illustrating data types.

  implicit none 
  real :: x
  real (kind=8) :: y, z
  integer :: m
  real (kind=8) :: i

  ! Comment out implicit none and see what happens
  ! with the following three lines when uncommenting them:
  !i = 10; print*, i
  !i = 10.2394872938479287; print *,i
  !a = 1.3234 ; print*,a
  
  m = 3
  print *, " "
  print *, "M = ",M   ! note that M==m  (case insensitive)


  print *, " "
  print *, "x is real (kind=4)"
  x = 1.e0 + 1.23456789e-6
  print *, "x = ", x


  print *, " "
  print *, "y is real (kind=8)"
  print *, "  but 1.e0 is real (kind=4):"
  y = 1.e0 + 1.23456789e-6
  print *, "y = ", y


  print *, " "
  print *, "z is real (kind=8)"
  print *, "  and 1.d0 is real (kind=8):"
  z = 1.d0 + 1.23456789d-6
  print *, "z = ", z

end program demo1

Download this code

Comments:

  • Exclamation points are used for comments

  • The implicit none statement in line 7 means that any variable to be used must be explicitly declared.

    Note

    Fortran has a very weird feature called implicit variables. With this feature, you don’t need to declare any variables! For instance, you can freely use any variables if you want, and variables name starts with i, j, k, l, m, n will be considered as an integer, and any other variables assumed to be a real number.

    This is an ancient feature for saving memory and punched cards. But we are living in the 21th century, so we don’t need to use this troublesome feature anymore!

    Thus, I highly recommend to you to use implicit none statement on every Fortran file and routine to turn off the implicit feature, and make all variables you are using as explicit. Or, you can use compiler flag like -fimplicit-none (in gfortran) to force the compiler use implicit none statement on every single routine.

  • Lines 8-10 declare four variables x, y, z, n. Note that x is declared to have type real which is a floating point number stored in 4 bytes, also known as single precision. This could have equivalently been written as

    real (kind=4) :: x

    y and z are floating point numbers stored in 8 bytes (corresponding to double precision in older versions of Fortran). This is generally what you want to use.

  • Fortran is not case-sensitive, so M and m refer to the same variable!!

  • 1.23456789e-10 specifies a 4-byte real number. The 8-byte equivalent is 1.23456789d-10, with a d instead of e. This is apparent from the output below.

Compiling and running this program produces

$ gfortran demo1.f90 -o demo1.exe
$ ./demo1.exe

 M =            3

 x is real (kind=4)
 x =    1.00000119

 y is real (kind=8)
   but 1.e0 is real (kind=4):
 y =    1.0000011920928955

 z is real (kind=8)
 z =    1.0000012345678899

For most of what we’ll do in this class, we will use real numbers with (kind=8). Be careful to specify constants using the d rather than e notation if you need to use scientific notation.

(But see Default 8-byte real numbers below for another approach.)

Intrinsic functions

There are a number of built-in functions that you can use in Fortran, for example the trig functions:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
! /codes/builtinfcns.f90

program builtinfcns

  implicit none
  real (kind=8) :: pi, x, y

  ! compute pi as arc-cosine of -1:
  pi = acos(-1.d0)  ! need -1.d0 for full precision!

  x = cos(pi)
  y = sqrt(exp(log(pi)))**2

  print *, "pi = ", pi
  print *, "x = ", x
  print *, "y = ", y

end program builtinfcns

Download this code

This produces

$ gfortran builtinfcns.f90
$ ./a.out
 pi =   3.1415926535897931
 x =   -1.0000000000000000
 y =    3.1415926535897927

See http://www.nsc.liu.se/~boein/f77to90/a5.html for a good list of other intrinsic functions.

Default 8-byte real numbers

Note that you can declare variables to be real without appending (kind=8) if you compile programs with the gfortran flag -fdefault-real-8, e.g. if we modify the program above to:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
! /codes/builtinfcns2.f90

program builtinfcns

  implicit none
  real :: pi, x, y    ! note kind is not specified

  ! compute pi as arc-cosine of -1:
  pi = acos(-1.0)

  x = cos(pi)
  y = sqrt(exp(log(pi)))**2

  print *, "pi = ", pi
  print *, "x = ", x
  print *, "y = ", y

end program builtinfcns

Download this code

Then

$ gfortran builtinfcns2.f90
$ ./a.out
 pi =   3.141593
 x =   -1.000000
 y =    3.141593

gives single precision results, but we can obtain double precisions with

$ gfortran -fdefault-real-8 builtinfcns2.f90
$ ./a.out
 pi =   3.1415926535897931
 x =   -1.0000000000000000
 y =    3.1415926535897927

Note that if you plan to do this you might want to define a Unix alias, e.g.

$ alias gfort8="gfortran -fdefault-real-8"

so you can just type

$ gfort8 builtinfcns2.f90
$ ./a.out
 pi =   3.1415926535897931
 x =   -1.0000000000000000
 y =    3.1415926535897927

Such an alias could be put in your .bashrc or .bash_profile (see Basic Unix/Linux Commands).

We’ll also see how to specify compiler flags easily in a Makefiles.

Fortran Arrays

Note that arrays are indexed starting at 1 by default, rather than 0 as in Python or C. Also note that components of an array are accessed using parentheses, not square brackets!

Arrays can be dimensioned and used as in the following example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
! /codes/array1

program array1

  ! demonstrate declaring and using arrays

  implicit none
  integer, parameter :: m = 3, n=2
  real (kind=8), dimension(m,n) :: A 
  real (kind=8), dimension(m) :: b 
  real (kind=8), dimension(n) :: x 
  integer :: i,j

  ! initialize matrix A and vector x:
  do j=1,n
     do i=1,m
        A(i,j) = i+j
     end do
     x(j) = 1.
  end do

  ! multiply A*x to get b:
  do i=1,m
     b(i) = 0.
     do j=1,n
        b(i) = b(i) + A(i,j)*x(j)
     end do
  end do

  print *, "A = "
  do i=1,m
     print *, A(i,:)   ! i'th row of A
  end do
  print "(2d16.6)", ((A(i,j), j=1,2), i=1,3)
  print *, "x = "
  print "(d12.4)", x
  print *, "b = "
  print "(d16.6)", b

end program array1

Download this code

Compiling and running this code gives the output

A =
  2.0000000000000000        3.0000000000000000
  3.0000000000000000        4.0000000000000000
  4.0000000000000000        5.0000000000000000
   0.200000D+01    0.300000D+01
   0.300000D+01    0.400000D+01
   0.400000D+01    0.500000D+01
x =
 0.1000D+01
 0.1000D+01
b =
   0.500000D+01
   0.700000D+01
   0.900000D+01

Comments:

  • In printing A we have used a slice operation: A(i,:) refers to the i’th row of A. In Fortran 90 there are many other array operations that can be done more easily than we have done in the loops above. We will investigate this further later.

  • Here we set the values of m,n as integer parameters before declaring the arrays A,x,b. Being parameters means we can not change their values later in the program.

  • It is possible to declare arrays and determine their size later, using allocatable arrays, which we will also see later.

There are many array operations you can do, for example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
! /codes/vectorops.f90

program vectorops

  implicit none
  real(kind=8), dimension(3) :: x, y

  x = (/10.,20.,30./)
  y = (/100.,400.,900./)

  print *, "x = "
  print *, x

  print *, "x**2 + y = "
  print *, x**2 + y

  print *, "x*y = "
  print *, x*y

  print *, "sqrt(y) = "
  print *, sqrt(y)

  print *, "dot_product(x,y) = "
  print *, dot_product(x,y)


end program vectorops

Download this code

produces

x =
  10.000000000000000        20.000000000000000        30.000000000000000
x**2 + y =
  200.00000000000000        800.00000000000000        1800.0000000000000
x*y =
  1000.0000000000000        8000.0000000000000        27000.000000000000
sqrt(y) =
  10.000000000000000        20.000000000000000        30.000000000000000
dot_product(x,y) =
  36000.000000000000

Note that addition, multiplication, exponentiation, and intrinsic functions such as sqrt all apply component-wise.

Multidimensional arrays can be manipulated in similar manner. The product of two arrays x and y using * (i.e., x*y) is always component-wise. For matrix multiplication, use matmul. There is also a transpose function:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
! /codes/arrayops.f90

program arrayops

  implicit none
  real(kind=8), dimension(3,2) :: a
  real(kind=8), dimension(2,3) :: b
  real(kind=8), dimension(3,3) :: c
  real(kind=8), dimension(2) :: x
  real(kind=8), dimension(3) :: y
  integer i

  a = reshape((/1,2,3,4,5,6/), (/3,2/))

  print *, "a = "
  do i=1,3
     print *, a(i,:)   ! i'th row
  end do

  b = transpose(a)

  print *, "b = "
  do i=1,2
     print *, b(i,:)   ! i'th row
  end do

  c = matmul(a,b)
  print *, "c = "
  do i=1,3
     print *, c(i,:)   ! i'th row
  end do

  x = (/5,6/)
  y = matmul(a,x)
  print *, "x = ",x
  print *, "y = ",y

end program arrayops

Download this code

produces

a =
  1.0000000000000000        4.0000000000000000
  2.0000000000000000        5.0000000000000000
  3.0000000000000000        6.0000000000000000
b =
  1.0000000000000000        2.0000000000000000        3.0000000000000000
  4.0000000000000000        5.0000000000000000        6.0000000000000000
c =
  17.000000000000000        22.000000000000000        27.000000000000000
  22.000000000000000        29.000000000000000        36.000000000000000
  27.000000000000000        36.000000000000000        45.000000000000000
x =    5.0000000000000000        6.0000000000000000
y =    29.000000000000000        40.000000000000000        51.000000000000000

Loops

Consider a code with do-loops:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
! /codes/loops1.f90

program loops1

  implicit none
  integer :: i

  do i=1,3           ! prints 1,2,3
     print *, i
  end do

  do i=5,11,2        ! prints 5,7,9,11
     print *, i
  end do

  do i=6,2,-1        ! prints 6,5,4,3,2
     print *, i
  end do

  i = 0
  do while (i < 5)   ! prints 0,1,2,3,4
     print *, i
     i = i+1
  end do

end program loops1

Download this code

The while statement used in the last example is considered obsolete. It is better to use a do loop with an exit statement if a condition is satisfied. The last loop could be rewritten as

i = 0
do                 ! prints 0,1,2,3,4
   if (i>=5) exit
   print *, i
   i = i+1
end do

This form of the do is valid but is generally not a good idea. Like the while loop, this has the danger that a bug in the code may cause it to loop forever (e.g. if you typed i = i-1 instead of i = i+1).

A better approach for loops of this form is to limit the number of iterations to some maximum value (chosen to be ample for your application), and then print a warning message, or take more drastic action, if this is exceeded, e.g.:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
! /codes/loops2.f90

program loops2

  implicit none
  integer :: i,j,jmax

  i = 0
  jmax = 100
  do j=1,jmax        ! prints 0,1,2,3,4
     if (i>=5) exit
     print *, i
     i = i+1
  end do

  if (j==jmax+1) then
     print *, "Warning: jmax iterations reached."
  end if

end program loops2

Download this code

Note: j is incremented before comparing to jmax.

if-then-else

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
! /codes/ifelse1.f90

program ifelse1

  implicit none
  real(kind=8) :: x
  integer :: i

  i = 3
  if (i<2) then
     print *, "i is less than 2"
  else
     print *, "i is not less than 2"
  end if

  if (i<=2) then
     print *, "i is less or equal to 2"
  else if (i/=5) then
     print *, "i is greater than 2 but not equal to 5"
  else 
     print *, "i is equal to 5"
  end if

end program ifelse1

Download this code

Comments:

  • The else clause is optional

  • You can have optional else if clauses

There is also a one-line form of an if statement that was seen in a previous example on this page

if (i>=5) exit

This is equivalent to

if (i>=5) then
   exit
end if

Booleans

  • Compare with <, >, <=, >=, ==, /=. You can also use the older Fortran 77 style: .lt., .gt., .le., .ge., .eq., .neq..

  • Combine with .and. and .or.

For example

((x>=1.0) .and. (x<=2.0)) .or. (x>5)

A boolean variable is declared with type logical in Fortran, as for example in the following code:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
! /codes//boolean1.f90

program boolean1

  implicit none
  integer :: i,k
  logical :: ever_zero
  real :: myTrulyTrulyVeryLongVariableNameToStoreRealVariable
  real :: a, b
  
  ever_zero = .false.
  do i=1,10
     k = 3*i - 1
     print*, i, k, ever_zero, (k == 0)
     ever_zero = (ever_zero .or. (k == 0))
  end do

  if (ever_zero) then
     print *, "3*i - 1 takes the value 0 for some i"
  else
     print *, "3*i - 1 is never 0 for i tested"
  end if

  a = 1.0
  myTrulyTrulyVeryLongVariableNameToStoreRealVariable = 2.0
  b = a + myTrulyTrulyVeryLong&
       &VariableNameToStoreRealVariable

  print*,a,b
end program boolean1

Download this code

Line Continuation

In case you wish to write a long line so that you want to continue your implementation in the next line, it can be continued using &

if (i>=5) &
   print *, 'This line is too long and I am using & to continue. This is handy!!!'

There is also a case you wish to end the line with & and begin a new line with &. For instance, if you have a very long variable name and want to split it off and continue in the next line (yes, this is a very weird coding practice though), you can do something like this

real :: myTrulyTrulyVeryLongVariableNameToStoreRealVariable
real :: a, b

a = 1.0
myTrulyTrulyVeryLongVariableNameToStoreRealVariable = 2.0
b = a + myTrulyTrulyVeryLongVariableNameToStoreRealVariable

or, the last can be put into two lines using &

b = a + myTrulyTrulyVeryLong&
        &VariableNameToStoreRealVariable

But, in this case, make sure that there is no space between the last chracter Long and the & that follows it

b = a + myTrulyTrulyVeryLong      &
        &VariableNameToStoreRealVariable

This will be then equivalent to

b = a + myTrulyTrulyVeryLong      VariableNameToStoreRealVariable

which is not what you wanted to do.

The situation is the same with the second line that follows the first line

b = a + myTrulyTrulyVeryLong&
    &      VariableNameToStoreRealVariable