Your Federal Quarterly Tax Payments are due April 15th Get Help Now >>

Tutorial introduction to Fortran 90 - PDF by msh36617


									    Tutorial introduction to Fortran 90

    A simple program

    The study of physics usually begins with mechanics, so it makes sense to start our study
    of computational physics in the same place.

    The intent of this tutorial is to allow you to get some familiarity with the programming
    process — creating a code file, compiling and debugging a program — and not to teach
    you to program. The program examples can be typed in as is, and along the way you will
    get some idea of the language.

    Of course, there are other good reasons. Mechanics has evolved methods of extreme
    sophistication in the 350 years since Newton. One of the main impetuses has been the
    study of orbital motion in the solar system. But these sophisticated methods are there to
    provide a good starting point for perturbation treatments; often obscured in the beauty
    of the formalism is the fact that there are few exact solutions. So numerical treatments
    — starting already with Euler — have provided a powerful alternative.

    Let’s start off with something very simple: dropping a ball from, say, the leaning tower
    of Pisa. The equation of motion is
                                               d2 y
                                           m        = −mg
    with y the height, t the time, m the mass of the ball and g the acceleration due to gravity.

    Of course, this has the solution at time t, assuming an initial height y0 and velocity zero
    at time t = 0,
                                          y = y 0 − g t2

    We can begin our programming in the simplest way, using the computer to evaluate this
    mathematical expression. We will evaluate the height for a set of times 0.1, 0.2, . . . , 1.0
    s. We will do this pretty much in the same way we would do the task with a scientific
    calculator, so think of the program as a list of instructions for using such a tool.

    Here is a simple program in Fortran to calculate the height as a function of time:

    program drop0
2   ! drop0.f90 - rwf
    ! Simple Fortran program
4   ! straightforward evaluation of expressions

6   implicit none

 8   real : :   g =9.8         !    accel due to gravity (m/s^2)
     real : :   m =1.0         !    mass of ball ( kg)
10   real : :   y0 = 1 0 . 0   !    initial height 10 m
     real : :   v0 = 0 . 0     !    initial vertical velocity component (m/s)
12   real : :   tmax = 1 . 0   !    drop time (s)
     real : :   dt = 0 . 1     !    time step
14   integer    : : Nsteps , step
     real : :   y
     ! determine number of time steps:
18   Nsteps=tmax / dt+1
     ! calculation loop
20   do step =1, Nsteps
         t = step ∗ dt
22       y = y0 −0.5∗ g∗t ∗∗2
         print ∗ , t , y
24   end do

26   end program drop0

     (Note that the line numbers printed here are just for reference below, and are not part of
     the program!)

        • Comments, which are there to provide explanations of what is going on, start from
          an exclamation mark and continue to the end of the line.
        • Note that the program is “bracketed” by program ... and end program ... state-
          ments. The name given in the statement is arbitrary, and does not have to be the
          same as the file.
        • implicit none is a statement that turns off some features of older versions of
        • Note that all variables used in the program are declared in the first section of the
          program; some are defined to have initial values.
        • Lines 8-15 define various constants and variables. Numerical values could simply
          be used instead in the appropriate places, but it is usually better to give names to
          these and not have “magic numbers” floating about in your program.
        • This program uses a do loop to iterate over a range of integers which are used to
          calculate the time and hence the height of the ball.
        • In line 18, the number of steps to be taken is calculated.
        • Line 20 introduces the do loop: The variable step is given values in succession from
          1 to Nsteps and the indented portion below is executed. Note that the statement
          ends with a statement end do. Do loops can be contained in other do loops, but
          cannot overlap.

   • Lines 21-23 constitute the body of the loop. Note that the indentation is not required,
     but is useful for seeing the structure of the code.

   • Note the exponentiation operator is **.

   • The print statement will print the values of the variables in its argument list on the
     terminal. The asterisk * is a placeholder for a format statement, and implies a free
     format. You can specify the format of the data by using a proper format statement
     instead of the asterisk; this allows you to specify number of significant digits, etc.
     Note the comma after the *.

Type this into a file. The comments are unimportant from the compiler’s point of view
and can be omitted — they are there for you to keep track of things. Save the file with a
.f90 extension, e.g. something like drop0.f90.

Now, before we can run the program, it must be compiled. We will use the program ifc
(Intel Fortran compiler) for this. The command for this is

ifc -o drop0 drop0.f90

This takes the file drop0.f90 and creates an executable program drop0 (that’s the -o

If the command prompt returns after this with no messages all went well. On the other
hand it is probable that you got some stream of error messages and warnings. You will
have to see what they say and try and fix your code, maybe with some help.

OK, it compiles cleanly. Now run it:


This line runs the program drop0 that exists in the current directory (that’s what the ./
indicates). (Technical note: the Unix PATH, which lists where executable programs exist,
usually excludes the current directory for security reasons, so you need to explicitly set
the path).

You program runs! You get a list of numbers: do they seem sensible?

Now make the problem a little less trivial by including the drag force, so that the equation
of motion is
                                  d2 y             1
                                m 2 = −mg + CρAv 2
                                   dt              2
Here, C is a constant; ρ the density of air and A the cross-sectional area of the body.

     This actually has the solution,

                                                      1            1 + (γv0 )2
                                         y = y0 +        ln
                                                     γ2g           cosh γgt

     where γ =        CρAv 2 /2mg and y0 is the initial height.

     This expression has a number of mathematical functions in it, including the square root.
     There are suitable versions of these supplied by the system. Here is a simple program in
     Fortran to evaluate this expression. You can enter this, but it’s quicker to copy your first
     program to a new file (cp drop0.f90 drop1.f90) and modify the new file. Use what
     you have.

     program drop1
 2   ! drop1.f90 - rwf
     ! Simple Fortran program
 4   ! straightforward evaluation of expressions

 6   implicit none

 8   real : :    g =9.8         !        accel due to gravity (m/s^2)
     real : :    m =1.0         !        mass of ball ( kg)
10   real : :    y0 = 1 0 . 0   !        initial height 10 m
     real : :    v0 = 0 . 0     !        initial vertical velocity component (m/s)
12   real : :    tmax = 2 . 0   !        drop time (s)
     real : :    dt = 0 . 1     !        time step
14   real : :    C =0.5         !        drag coefficient
     real : :    rho = 1 . 0    !        density of air ( kg/m^3)
16   real : :    R =0.1         !        radius of ball (m)
     real : :    pi =3.1415926
18   integer     : : Nsteps , step
     real : :    y , t , gamma
     ! determine number of time steps:
22   Nsteps=tmax / dt+1
     ! calculation loop
24   do step =1, Nsteps
         t = step ∗ dt
26       gamma=sqrt ( C∗ rho ∗ pi∗R ∗ ∗ 2 / ( 2 ∗ m∗g ) )
         y = y0 + ( 1 . 0 / ( g∗ gamma ∗ ∗ 2 ) ) ∗ log ( sqrt ( 1 . 0 + ( gamma ∗v0 ) ∗ ∗ 2 ) / ( cosh ( gamma ∗g∗t ) )
28       print ∗ , t , y
     end do
     end program drop1

         • In many languages you will need to import the math functions from a library or
           module. Fortran has them built into the language. (You can of course import other
           things from a library or module).

       • Note the way of calling a function — this is very similar to obvious mathematical
         notation. A function typically returns a value.

       • Note that trigonometric functions expect their argument in radians.

       • The system does not supply a value for π.

    The distance fallen in a given time is presumably smaller than in the original code, where
    there was no air resistance. How much smaller? Calculate it! Let’s add a line to calculate
    the vacuum value of free fall, say yvac = y0-0.5*g**2. Add this line to the do loop, and
    alter your program to print it out. Note that you have a new variable name.

    How about a graph? The simplest way to do this is to use a program like gnuplot. This
    requires the data to be plotted to be in a file. The simplest way to do this is to cut and
    paste from the screen to a text editor, and save the file, for example as drop1.out. The
    start gnuplot and plot it:

    gnuplot> plot ’drop1.out’

    Hmm. This only plots one set of data. You have to tell the program which columns to

    gnuplot> plot ’drop1.out’ using 1:2, ’drop1.out’ using 1:3

    This can be abbreviated; you might like to use lines instead of points.

    gnuplot> plot ’drop1.out’ u 1:2 with lines, ’drop1.out’ u 1:3 w l

    Now let’s change the program so that it runs until one of the objects hits the ground at
    y = 0. We can do this by replacing the do loop with a do-while loop:

    do while ( yvac > 0 . 0 )
2       step = step + 1

    While we’re at it change the distance fallen to 50 m.
    Now, when we run the program, the data scrolls out of view. It’s time to write it to a file
    from within the program. We can do so by first opening a file

    open(1, file=’drop2.out’)

Note that the file name is a character string enclosed in quotes. the “1” is a logical unit
used to refer to the file. To write to the file we can replace the print statements with:

write(1,*) t,y,yvac

Here the “1” stands for the file; the “*” is still a placeholder for a format. Note that there
is no comma before the t. Finally, we can close the file (after we have written everything
to it.


Plot the results using gnuplot.


To top