MODULE 4A2 INTRODUCTION TO FORTRAN
FORTRAN was the first high level programming language to come into general use and it is still
by far the most widely used language in Mechanical, Civil and Aeronautical Engineering. The
advantage of FORTRAN is that it gives very fast executable code, which minimises the run time
for long "number crunching" calculations.
The language has evolved considerably over the years and the latest version, FORTRAN95,
includes many of the features of more modern languages such as "C" . However, for our purposes
these features merely make things more complicated and we will use an older version
FORTRAN77 which is still very widely used.
By far the best way to learn the language is to study existing programs and to start to write one
yourself. The following notes are intended to give only a very brief overview of the language.
Any textbook on FORTRAN77 will give far more details.
A FORTRAN program consists of a main program and any number of subroutines. The
Subroutines are really separate programmes that are called by the main program to perform
specialist operations. Vast libraries of subroutines are available to perform such operations as
Matrix Inversion, Interpolation, Fourier Transforms, Statistical Analysis, etc. However, in many
cases the subroutines are written by the developer of the program and are used as a way of
simplifying the layout of the main program. The layout of a typical program might be:
C TITLE OF PROGRAM
The line beginning with C in the first column are just comments and are not executed. Executable
statements start in column 7.
A subroutine is called from the main program by a call statement such as
And the subroutine layout is
Data is transferred from the main program to the subroutines and back again by either "Common
Blocks" or by calling the subroutine with arguments. In the former case both main program and
the subroutine contain statements like:
COMMON X,Y,Z, TIME, A(100),----
This means that the variables X,Y,Z,TIME, ---etc, are available for use in both the main program
and the subroutine.
In the latter case the call from the main program to the subroutine would be something like
CALL SUB1(X,Y,Z,TIME,A ---- )
And the subroutine would be headed
SUBROUTINE SUB1(X,Y,X,TIME,A---- )
The variables are then transferred from the main program to the subroutine and back to the main
program when the subroutine has completed its operations.
Variables may be either simple variables, X,Y,Z, I,J,K, etc . Or they may be arrays like
X(100),Y(150),A(100,100), K(1000), ---etc. The latter means that there are 100 values of X,
X(1),X(2),X(3), ---- X(100), etc .
Variables may be used without declaring their type, real or integer, because BY DEFAULT any
variables beginning with I,J,K,L,M,N are taken to be integers, all other variables are taken to
be real. The default can be overridden by formal declaration of the type of a variable, e.g.
Would make MACH a real number and STEP an integer. However, in most cases it is more
convenient to use the defaults.
A very common source of errors is to forget this convention and use a variable starting with
I,J,K,L,M,N as a real number of vice-versa.
Any variables, real or integer, which are to be used as arrays must be declared by a dimension
statement giving the size of the array. This must be done at the start of the program before the
array is used. e.g.
DIMENSION XX(100),YY(200), AREA(100,100), NSTEP(1000), --etc.
If an array is transferred as an argument of a subroutine then its dimensions must also be declared
in the subroutine, see the subroutine call example in the last section.
The array dimensions may also be declared in a COMMON statement,
COMMON XX(100),YY(200), AREA(100,200), NSTEP(1000), --etc.
Which must appear in any subroutine where any of the variables in the common block are used.
Variables cannot be transferred from the main program to a subroutine by both a common block
and as subroutine arguments. Only one or other method may be used.
LAYOUT OF STATEMENTS
FORTRAN programmes were first written on punched cards and so the layout of statements is a
legacy from that.
The statements may use columns 1-72 of the page or screen. However, columns 1-5 are reserved
for a label, i.e. a number referencing the statement.
Any character in column 6 means that the line is a continuation of the previous statement.
The executable statements proper only start on column 7 and may continue to column 72. A
"C" in column 1 means that the line is merely a comment to help understanding of the program
and is not a FORTRAN statement. Such comments greatly help to improve program layout and
understanding. Blank spaces on a line are ignored and may again be used to improve layout.
Any characters beyond column 72 are ignored. This is a very common source of error.
FORTRAN does not distinguish between upper and lower case characters. Preferences vary so
choose your own. However, some computers do distinguish upper and lower case in file names.
File "DATA1" may not the same as file "data1".
C EXAMPLE OF SIMPLE FORTRAN PROGRAM TO CALCULATE THE
C AREA OF A CIRCLE
OPEN(UNIT = 6, FILE='RESULTS')
PI = 3.14159
DO 100 J=1,100
R(J) = J
AREA(J) = PI*R(J)*R(J)
WRITE(6,*) 'PROGRAM COMPLETED'
All arithmetic operators must be used explicitly. e.g. X = (A+B)*C
NOT X = (A+B)C
+ , - are as in ordinary arithmetic.
* is used or multiplication and / for division.
** is used for exponentiation , e.g. x**3 is x cubed. Y**2.5 is Y to the power 2.5.
Operations in brackets are performed first, followed by multiplication and division in the order
they appear on the line and by finally addition and subtraction
A/B*C means (A/B)*C NOT A/(B*C)
A + B/C is A + (B/C) NOT (A+B)/C
Arithmetic operations may mix real and integer numbers, e.g. Y = X*K . However, the result is
integer or real as determined by the type of the variable name, e.g. J = X*K will make J an
integer which will be truncated from the exact value of X*K whilst Y = X*K will make Y a
real number, the exact value of X*K.
A very common source of error is to divide two integers and set the result to a real number,
e.g. X = J/K . The value of X will be truncated to the nearest integer below the exact result of
dividing J by K. e.g. If J = 2 and K=3 the above will give X= 0. To obtain the correct result you
must first convert the integers to real numbers and then do the division i.e.
X = J/K gives X = 0.
X = FLOAT(J)/FLOAT(K) gives X = 0.6666666 .
This is a very common source of error and can waste hours of your time.
There are library functions available for most common algebraic functions. The argument of the
operator is always placed in brackets. e.g. SIN(X), COS(Y), TAN(W) for which the argument
must be in radians.
Other functions are: SQRT, LOG ( = loge), LOG10 , ATAN, ASIN, ACOS , EXP , ABS,
SINH, COSH, TANH , etc. A complete list can be found in any FORTRAN book.
Note that the ASIN, ACOS and ATAN functions always return values between 0 and Pi. This
can cause problems when a result between Pi and 2*Pi is required. This can be overcome for the
ATAN function by using ATAN2(Y,X) which, given values of Y and X, always gives the true
angle between the origin and the point with coordinates (X,Y). i.e. ATAN2(y,x) = arctan(y/x).
e.g. ATAN2(1.732,1) = pi/3 . This may be useful in your CFD exercise.
It is very common to perform the same operation many times using different values of the
variables. This is best done by a loop such as:
DO 100 J=1,JMAX
A(J) = PI*R(J)**2
which calculates JMAX values of A(J) . This can equally well be written :
A(J) = PI*R(J)**2
For long loops it is probably better to use the former layout since it is easier to locate the end of
Multiple nested loops are common in CFD calculations, e.g.
DO 200 K=1,KMAX
DO 200 J=1,JMAX
DO 200 N=1,NMAX
A(N,J,K) = B(N,J,K)*C(K,J)*R(N)
These use the logical operators GT, EQ, LT, GE, LE, NE, AND, OR to check whether a
statement is true or not and switch the subsequent operations accordingly. Usually the layout is
IF(A.GT.B) GO TO 100
or IF(A.GT.B.AND.N.LT.NMAX) THEN
or IF(EROR.GT.EMAX) THEN
EMAX = EROR
EMAX = 0.0
INPUT AND OUTPUT FILES
Input and output to a FORTRAN program is usually read from files or to/from the screen. Each
file that is to be used must be given a number, which is called its UNIT, and a name. The name of
the file is the same as the name that you give it when editing it. The file (UNIT) must be opened
before the file can be read or written to. This is done as follows:
An input file from which data is being read must exist and contain the required data before it is
opened. An output file need not exist before output is written to it, it will be created by the
FORTRAN program when it is run.
A special name is reserved for the screen , '/dev/tty'
To send output to the screen and call it unit 5 you would open it by
Once you have finished reading or writing a file it is safer to close it. To close UNIT 1 you write.
Once closed the file cannot be used again until it is reopened.
By default UNIT 6 is taken to be the screen and so you can write output to UNIT 6 without
bothering to open it. This is often useful for temporarily writing out variables whist debugging a
READING AND WRITING INPUT AND OUTPUT
Reading and writing input and output neatly using FORMAT statements can be the most difficult
part of FORTRAN. We will avoid most of the difficulties by always using unformatted input and
output. This may not look tidy but it is perfectly adequate when the results are being presented
graphically, as we will usually do.
For unformatted input and output we use: READ(N,*) and WRITE(N,*) where N is the number
of the UNIT (file) we wish to read from or write to. The numbers in the input file are read in
sequence with one or more spaces between adjacent numbers. Similarly the numbers in the
output file are written in sequence but the spacing may be erratic and the output may not be very
e.g. READ(1,*) X,Y,Z
To read or write text it is most easy to simply place it in single quotation marks and then treat it
exactly as numerical data, e.g.
WRITE(1,*) ' INPUT THE NUMBER OF POINTS IN THE DATA SET '
WRITE(1,*) ' STOPPING BECAUSE NPOINTS > 100 .'
READ(2,*) (R(N), N=1,NPOINTS)
DO 100 N = 1,NPOINTS
AREA = PI*R(N)*R(N)
CIRC = 2.*PI*R(N)
WRITE(1,*) ' PROGRAM COMPLETED '
Writing graphical output is much more complex and depends on the graphics package being
used. In effect a graphics package contains a library of subroutines that you can call from your
program. In 4A2 we use a package called HGRAPH written by Prof. Hodson at the Whittle lab.
However, graphics programs have been written to present all your results and there is no need for
you to learn how to use this package.
THE "INCLUDE" STATEMENT
Often we want to include the same piece of code in several different subroutines. To avoid
having to type it in every time it is needed it is possible to write it as a separate file and use the
INCLUDE statement to insert it wherever it is needed. This is particularly useful for inserting
common blocks because it avoids errors which will occur if the blocks are in any way different
in the different subroutines.
To do this we call the lines of code to be inserted by any file name, e.g. INSERT and whenever
we wish to include it in out program we use:
Where the file name is enclosed in single quotes.
e.g. if the file INSERT is the two lines
& NMAX,IMAX,JMAX, ALPHA
We can include it in SUBROUTINE UPDATE by using:
DO 100 N=1,NMAX
VX(N) = VY(N)*COS(ALPHA)
DO 50 I=1,IMAX
DO 50 J=1,JMAX
AREA(I,J) = PI*R(I,J)**2
Note how the variables NMAX, VX, VY, ALPHA , etc are immediately available for use in the
It will be exceptional if during the 4A2 exercises you do not spend several fruitless hours trying
to find BUGS in your code. The FORTRAN compiler in the DPO is not very helpful and you
have to locate the errors by studying your code in great detail and writing out intermediate
It greatly helps debugging if your code is laid out neatly and if long
expressions are broken down into several simpler expressions.
e.g. to find the velocity components induced by a vortex strength GAMMA located at X1, Y1 at
the point X2, Y2 we could write
READ(1,*) GAMMA, X1,Y1,X2,Y2
VX = GAMMA/(2*3.14159*SQRT((X2-X1)**2 + (Y2-Y1)**2))*(Y2-Y1)
& /SQRT((X2-X1)**2 + (Y2-Y1)**2)
VY = -GAMMA/(2*3.14159*SQRT((X2-X1)**2 + (Y2-Y1)**2))*(X2-X1)
& /SQRT((X2-X1)**2 + (Y2-Y1)**2)
(Try checking if the number of brackets is correct !! )
But it is much clearer and much easier to check if you write:
READ(1,*) GAMMA, X1,Y1,X2,Y2
DX = X2-X1
DY = Y2-Y1
DRSQ = DX*DX + DY*DY
2PI = 2.0*3.14159
VX = GAMMA*DY/(2PI*DRSQ)
VY = -GAMMA*DX/(2PI*DRSQ)
To help debugging it is also very important to use variables whose names are easily
recognisable by both you and by the demonstrators. e.g. call velocity components VX and VY
or U and V not ALPHA and OMEGA.
It will greatly help both you and the demonstrators if you try to adopt this type of programming
COMMON PROGRAMMING ERRORS
There are many possible errors in coding, some of them are extremely subtle and can take an
experienced demonstrator hours to find. However, probably the most common errors are listed
below. Always look through this list to see if you can identify your error before calling a
1) Using an integer (I,J,K,L,M,N) when you want a real number. The variable will be
rounded down to the nearest integer. e.g. KE = 0.5*V**2.
2) Using an integer divide when you want the answer to be a real number. e.g. X = J/K .
This is the most common mistake made by students in the 4A2 module.
3) Using a symbol that you have declared to be an array as a single variable. HO = Cp*T,
where HO has been declared as HO(100,100).
If your program runs very slowly this is probably what you have done because the
compiler will make the program perform the operation for all the 10000 Values of HO,
not just the one you intend.
4) Using ACOS,ASIN or ATAN to obtain an angle, then taking a different trig function of
the angle. The angle found may be in the wrong quadrant so the second trig function may
have the wrong sign. Try to avoid ASIN, ACOS, ATAN and use ATAN2(y,x) if
5) Using a variable whose value has not been previously set either by reading it in as data or
by putting it on the LHS of an expression. Make sure that you read in or calculate all
variables before using them.
6) Confusing I and 1, especially in arrays A(1,J) may be typed instead of A(I,J) .
7) Continuing a statement beyond column 72. e.g. if a line ends in the variable AREA and
the first A falls in column 72 then the REA will be ignored and the variable A will be
used instead of AREA. Even if the variable A has not been set and is not used elsewhere
in the program the compiler will not complain and the value of A will be taken as zero.
COPYING, COMPILING AND RUNNING PROGRAMS
Three short programs are provided for you to study in order to help you to become familiar with
FORTRAN. Simple modifications to the programs are suggested for you to try yourself in order
to gain experience in FORTRAN programming. Listings of the programs are provided at the end
of this handout.
These programs are intended only as an aid to learning FORTRAN. You will not be marked on
your attempts to modify them and you should not spend too much time on it. In particular you
should not delay starting on the panel method exercise beyond one week after the first lecture.
To copy the programs and data follow the instructions below exactly, using lower case for all
Logon in the normal way.
Type start to start the windows manager.
Open a new window.
Create a new directory for your work on this module by typing:
mkdir c2dir or you can use any other name you wish instead of c2dir.
Move to this directory by typing cd c2dir.
Now copy the programs and data by typing:
cp /public/teach/4A2/projectile.f .
cp /public/teach/4A2/data1 .
cp /public/teach/4A2/comp-flow-tables.f .
cp /public/teach/4A2/vortex.f .
cp /public/teach/4A2/vortex.dat .
cp /public/teach/4A2/compile .
The dots at the end of each command are essential.
The command compile has been created to enable you to compile the FORTRAN programs in
order to produce executable code. It automatically links in the graphics library for future use.
Having copied the programs you can now compile them by typing, for the first program:
The executable code produced by compile is always named go and so to run the program after
compiling it you simply type: go
If you compiled projectile.f this should produce a picture on the screen showing the path of
When you compile another program then go will be overwritten by the new executable. If you
wish to save the old one you can rename it by typing
mv go projectile.x or any other recognisable name.
SUGGESTED MODIFICATIONS TO THE PROGRAMS
PROGRAM projectile.f .
This calculates the motion of a projectile fired at a specified velocity and angle in the
atmosphere. The program calculates the acceleration of the projectile, including the drag, and
plots the trajectory on the screen. It reads a data set data1 giving the initial velocity, angle, etc.
At present the projectile stops when it hits a horizontal surface through the launching point. It is
suggested that you modify it to make it stop when it hits a surface through the launching point but
sloping at an angle to the horizontal. Modify the data input to read in the angle of the slope and
modify the program to make the calculation stop when the projectile hits this surface.
PROGRAM comp-flow-tables.f .
This calculates some of the functions of Mach number which are tabulated in the compressible
flow tables and writes them to a file called tables . It is suggested that you modify the program
to evaluate the function of Mach number: 0.5 ρ V2/(Po-P) which gives an estimate of the
error in using Bernoulli's equation. You already have V/√Cp To, P/Po and T/To available in the
program so it is not difficult to calculate the above function.
Evaluate this within the "do loop". To print it out neatly using the FORMAT statements it is
easiest if you replace one of the other functions that are already being printed by the new
This program is of some relevance to the exercises that will follow. It calculates the motion of a
cloud of vortices each moving under the influence of the others. Remember that a vortex moves
with the fluid adjacent to it and does not induce any velocity on itself. The number of vortices,
number of steps, vortex strengths and positions are read in from a file vortex.dat and the paths of
all the vortices are plotted on the screen.
To modify the program it is suggested that you try to add the velocity due to a single sink, of
strength SINK fixed at the origin, i.e. the source does not move. Read in the value of SINK
from the file vortex.dat and add the velocity it induces at every vortex to the previous value of