In the first year you had your first taste of computing and programming. This workshop is designed to give you more practice and experience in using Fortran 90. Before you start please read this script carefully. It is important that you attempt the programming individually, not in pairs or as a consortium. You can only become proficient at programming by doing it yourself, not by watching someone else. The workshop is not a race and you must not worry if you seem to be falling behind others. The person who completes first will almost certainly have made mistakes or cut a few corners. You MUST keep a record of what you do in your laboratory workbook, obtain listings of all your programs and retain the program source files on your directory at the end of the 4 laboratory sessions so that your work can be accessed.
You should use the PCs situated between the 1st and 2nd Year Laboratories during timetabled sessions. Other facilities like the Bennett User Area and terminals at the Halls of Residence which provide connection to the Irix System can be used outside normal working hours. Hardcopy printouts of your source listings and graphs can be collected from the Print Station situated in the terminal room between the First and Second Year Laboratories. You will need a Print Copycard which can be obtained from the demonstrator.
There are many good books on Fortran 90. Notes on the syntax were provided last year and are extended below but they are not intended to be a comprehensive text in competition with published material. You don't have to buy a book but the following books are recommended. Whichever book you choose make sure it is up to date and covers the 1990 standard. Most books on Fortran 90 do contain a description of the 77 standard but be warned there are significant differences between f77 and f90. Do not purchase a book which only covers the 77 standard.
The script for this workshop, the scripts for other workshops, useful references on Fortran 90, documentation on the subroutine library PGPLOT and a guide to using Irix and Word can all be found on the World Wide Web at URL:
http://www.star.le.ac.uk/~rw/compshop/
You will need a user name and password for both the CFS system and the Irix service. First of all log on to Windows 2000. Then using the ``Start'' menu at the bottom left hand corner of the screen goto ``Programs'' and then ``Exceed'. If you start the Exceed program you will be given the choice of two machines, harrier and owl. Use the mouse to click on either harrier or owl. This will bring up a login screen for an Irix machine. Type in your user-name followed by <return>, then type your password and again complete the line by typing <return>. The <return> signals the end of a line to the machine and tells it to consider what you have typed. Note that UNIX (IRIX) is generally case sensitive and it is not advisable to use the caps lock key. You should now be logged-on and the prompt $ will appear on the left-hand side of your screen. This is generated by the operating system IRIX (a flavour of UNIX) soliciting for commands. When commands are described in the following text the computer types the prompt $ and you should type what follows terminated by <return>. To log off the machine type exit <return>
$ exit
After logging out of Irix you will drop back to the Windows 2000 machine. Don't forget to log off this using ``Logoff'' in the ``Start'' menu. Use ``Logout'' rather than ``Shutdown'' because the latter will require the machine to reboot which takes some time.
Please don't leave your terminal logged on and unattended.
This is fully explained in the first year workshop. If you have forgotten how to do this now is the time to look it up.
To summarise, use the following commands to create and run a program
$ eve <name>.f90 $ f90 <name>.f90 -o <name> $ <name>
You may prefer to use the nedit editor in place of eve. This is best run as a background task using the & qualifier
$ nedit <name>.f90 & $ f90 <name>.f90 -o <name> $ <name>If you want to find more about the UNIX operating system you can browse through the help manual using the man command.
$ man -k text $ man command
This uses the command more to give pages of information. The -k option allows you to search the manual for specific text sequences. Otherwise you must search for a specific command.
We are going to use a dynamics problem involving a magnetically perturbed pendulum as a vehicle to introduce the scope and power of F90. The problem cannot be solved analytically and exhibits complicated chaotic behaviour as well as simpler regular motion.
The magnetically perturbed pendulum in question is similar to some executive toys you can buy to adorn your desk and distract you when the going gets tough. It consists of a magnetic bob suspended on a wire or string which lies along the z-axis when the wire is vertical. The bob is free to move in both x and y and we will only consider the motion of the bob in the horizontal xy plane. To make the motion more interesting the bob is perturbed by a magnet which can be placed in the xy plane just underneath the swinging pendulum. When the bob passes near to the magnet the path is deflected violently and makes the pendulum swing in a random way. The geometry of the pendulum is illustrated in figure 1.
The motion is described by 2 simple coupled differential equations:
|
|
In each case the 3 terms on the right correspond to gravity, the magnetic perturbation and air drag. The constants g, m and d depend on the mass of the bob, the strength of the magnet and the viscocity of air. The fixed magnet is at position (xm ,ym , zm) and rm=((x-xm)2+(y-ym)2+(z-zm)2)1/2. Actually these equations are only an approximation to the real state of affairs but as you will be able to demonstrate for yourselves this doesn't matter.
If m and d are zero then the pendulum executes simple harmonic motion in x and y and the bob follows an elliptical path in the xy plane. If d is zero then the motion is conservative (no energy loss) and we can express the force using a potential energy:
|
where r=(x2 + y2)1/2. The forces in x and y (not including drag) are then given by the gradients:
|
|
The kinetic energy is given by:
|
The state of the bob at any time is specified by 4 numbers ( x, y, dx/dt , dy/dt). These define what is called a phase space for the system. As time progresses the bob will trace out a path in phase space. If we take a particular starting position for the bob in phase space we can calculate the path by integrating the differential equations. There we have it. The problem is to integrate the differential equations given g, m , d, and a starting position and to find the path of the bob in phase space.
In the context of Fortran programming we will represent the position in phase space by an array of 4 numbers, P(4) where P(1) and P(2) are the position in xy and P(3) and P(4) are the corresponding velocities.
In the following script the tasks you are expected to do are described in subsections entitled Task n.
Last year you learnt how to write a simple PROGRAM section. This year you will be exposed to the full power of F90; in particular we are going to use Modular Programming involving MODULEs, SUBROUTINEs and FUNCTIONs. Correct use of these makes the code easier to write, easier to maintain and alter and in many cases more efficient.
We will start with functions. These are just like the intrinsic functions SIN, COS and so forth but you define what is required and how to calculate the answers.
Suppose we have the position of the pendulum bob in phase space and want a function to return the kinetic energy. The code for the function might look like the following:
FUNCTION PEND_KE(P)
! Return the kinetic energy of pendulum
IMPLICIT NONE
INTEGER, PARAMETER :: DOUB_PREC=SELECTED_REAL_KIND(P=10,R=30)
REAL(DOUB_PREC) :: PEND_KE,P(4)
! Calculate kinetic energy using velocities
PEND_KE=0.5*(P(3)**2+P(4)**2)
END FUNCTION PEND_KE
The first line starts the definition of the function. It defines the name, in this case PEND_KE, and then defines the argument list, in this case just one variable P. The line starting INTEGER, PARAMETER defines a parameter DOUB_PREC which is used in the declarations of the REAL variables. We will say more of this later. Notice that the function name PEND_KE and the array P(4) are both declared REAL(DOUB_PREC). This is because the answer is returned in the name of the function as shown in the assignment line. Finally the line starting END completes the definition of the new function.
In order to use the function we must evoke it in a program:
PROGRAM TEST_FUN
IMPLICIT NONE
INTEGER, PARAMETER :: DOUB_PREC=SELECTED_REAL_KIND(P=10,R=30)
REAL(DOUB_PREC) :: PEND_KE,P(4),A
! Read in position in phase space
PRINT *,"type position x,y,vx,vy"
READ *,P
A=PEND_KE(P)
PRINT *,"Kinetic energy",A
END PROGRAM TEST_FUN
In the above you had to type the definition line for the DOUB_PREC parameter in both the function and the program module. It is often the case that the same declaration lines are required to appear in a large number of places. Furthermore it is often useful to define some variables as global between a number of routines rather than passing them as arguments. The MODULE provides a clean and uniform way to achieve these and other functions.
The following module defines various things which will be useful in the programming solution to our pendulum problem.
MODULE PEND_DEF
! Definitions for magnetic pendulum
IMPLICIT NONE
SAVE
! Select precision for real numbers
INTEGER, PARAMETER :: DOUB_PREC=SELECTED_REAL_KIND(P=10,R=30)
INTEGER, PARAMETER :: SING_PREC=SELECTED_REAL_KIND(P=6,R=30)
! Definition of data type for pendulum parameters
TYPE PEND_PARAM
! Position of fixed magnet
REAL(DOUB_PREC), DIMENSION(3) :: POS
! Force constants for magnet, gravity and drag
REAL(DOUB_PREC) :: FMAG,FGRA,FDRA
END TYPE PEND_PARAM
! Declare global for magnet parameters
TYPE(PEND_PARAM) :: G ! parameters of fixed magnet etc.
END MODULE PEND_DEF
The first line starts the definition and gives the module a name and the last line starting with END completes the definition. The IMPLICIT NONE and SAVE lines are recommended in all modules since they ensure that every variable is defined and is non- volatile. The INTEGER declarations define two integer constants which refer to REAL types of a specified precision (P=) and range (R=). These are required because single precision (P=6 decimal digits) is insufficient to give adequate precision for the solution to our problem. However when we come to plot the results later we must use single precision variables for the plotting routines. So all variables which are used in calculation will be declared REAL(DOUB_PREC) while plotting variables must be declared REAL(SING_PREC). The remaining lines in the MODULE PEND_DEF define and declare a new data type. In Fortran 90 new data types can be defined in term of existing data types. In this case a new type called PEND_PARAM consisting of a position for the fixed magnet (an array P(3)) and three constants FMAG, FGRA and FDRA (corresponding to m, g, d above). The lines between TYPE PEND_PARAM and END TYPE perform this definition. The new data type is then used to declare a (global) variable G which will hold all the required information about the pendulum.
That completes the definition of the module. To include the definitions provided by the module in a PROGRAM or FUNCTION (or SUBROUTINE, see later) you must put USE name as the first line after the start of the new routine.
As an example of this we are now going to produce a FUNCTION which returns the potential energy of the pendulum. If you refer back to the pendulum equations you will see this requires the constants m and g (FMAG and FGRA components of G). You can see below how the MODULE PEND_DEF is included into the FUNCTION PEND_PE and how the components of G are referred to using the syntax G%FMAG to pick out the value of m and G%POS(1) to pick out the x position of the magnet.
FUNCTION PEND_PE(P)
! Return the potential energy of pendulum
USE PEND_DEF
IMPLICIT NONE
REAL(DOUB_PREC) :: PEND_PE,P(4)
! Local variables
REAL(DOUB_PREC) :: R2,RM2,XM,YM,RMZ2
! Calculate radial coordinates
XM=P(1)-G%POS(1)
YM=P(2)-G%POS(2)
R2=P(1)**2+P(2)**2
RM2=XM**2+YM**2
RMZ2=RM2+G%POS(3)**2
! Calculate potential
PEND_PE=0.5*G%FGRA*R2-G%FMAG/(RMZ2**3)/6.0
END FUNCTION PEND_PE
PROGRAM TEST
USE PEND_DEF
IMPLICIT NONE
REAL(DOUB_PREC) :: PEND_KE,P(4),A
! Get magnet parameters
PRINT *, "x,y,z,mag,grav,drag"
READ *, G
...
Your source file should now contain the MODULE PEND_DEF, the FUNCTION PEND_KE, the FUNCTION PEND_PE and the PROGRAM TEST.
The program should read in the magnet parameters and the position of the bob in phase space and then test the two functions.
In order to solve the coupled differential equations of our pendulum problem we need to translate the equations into Fortran. This is best done using a SUBROUTINE which is the last form of module supported by F90. Unlike a FUNCTION all the information to and from a SUBROUTINE is passed as arguments. In general some arguments will be purely input, some will be input, modified by the routine and then output and others will be just output. We want a routine which will calculate the derivatives of the position in phase space given the position and time. The following subroutine does this:
SUBROUTINE PEND_EQM(T,P,PDOT)
! Differential equations of motion for the magnetic pendulum
USE PEND_DEF
IMPLICIT NONE
REAL(DOUB_PREC) :: T ! time
REAL(DOUB_PREC) :: P(4) ! position in phase space
REAL(DOUB_PREC) :: PDOT(4) ! returned derivatives
! Local variables
REAL(DOUB_PREC) :: XM,YM,RMZ,FM
! Calculate radius from fixed magnet
XM=P(1)-G%POS(1)
YM=P(2)-G%POS(2)
RMZ=SQRT(XM**2+YM**2+G%POS(3)**2)
! Calculate magnetic force
FM=G%FMAG/(RMZ**8)
! Calculate derivatives
PDOT(1)=P(3)
PDOT(2)=P(4)
PDOT(3)=-G%FGRA*P(1)-FM*XM-P(3)*G%FDRA
PDOT(4)=-G%FGRA*P(2)-FM*YM-P(4)*G%FDRA
END SUBROUTINE PEND_EQM
Read this carefully and convince yourself that it faithfully reproduces the required equations. The arguments T and P are the time and position (input arguments) and PDOT holds the output values of the derivatives. Like the functions we defined above the new routine uses the module PEND_DEF to provide definitions of DOUB_PREC and the global variable G.
REAL(DOUB_PREC) :: T ! time
REAL(DOUB_PREC) :: PDOT(4) ! returned derivatives
a print and read to get the time:
PRINT *,"Type time"
READ *,T
and the call to the subroutine followed by a print to list the result:
CALL PEND_EQM(T,P,PDOT)
PRINT *,"derivatives",PDOT
In order to perform the integration another subroutine is required. This routine starts as follows:
SUBROUTINE RUNKUT(DT,F,TIME,NV,P)
The routine performs an integration step by the fourth order Runge-Kutta method. You don't need to know what this is , you just need to know how to call the routine. The arguments of the routine are as follows:
REAL(DOUB_PREC) :: DT ! The time step EXTERNAL F ! The subroutine to return derivatives REAL(DOUB_PREC) :: TIME ! The time INTEGER :: NV ! The number of variables REAL(DOUB_PREC) :: P(NV) ! The new position in phase space
For the present application F will be the subroutine PEND_EQM and NV=4. RUNKUT will use subroutine PEND_EQM to predict the new position of the pendulum at time TIME+DT. The routine updates the TIME and P variables every call.
$ cp /disk/s/zrw/under/runkut.f90 . $ cat runkut.f90 your.f90 > your_new.f90
If you look at the routine you will see that variable names are not the same as above. This doesn't matter since variable names are only local to each routine unless they appear in a MODULE. In fact the RUNKUT subroutine has no connection with the pendulum problem other than it provides a method of solving ordinary differential equations. All the information about the problem is passed to the routine in the arguments F and NV.
The program required to integrate the pendulum equations must do the following:
In order to test the program is working you should calculate the answer in cases when you know the expected results. For example if d=0 and m=0 then the bob should execute simple harmonic motion. If g=1 what is the frequency? Choose the time increment DT such that you get several tens of samples for each complete cycle. Make sure that the motion gives the frequency as expected if g is altered.
DO J=1,NSTEPS
...
IF(MOD(J,100).EQ.0) THEN
... print stuff out every 100th pass
ENDIF
ENDDO
The best way of presenting many results is to plot a graph. If you plot x as a function of t as calculated above it should be easy to see if your results are sensible or not. If you plot x vs. y or x vs. vx then it is easier to visualize the path of the pendulum in phase space. In order to do this you need to use some form of graphics software. This is provided in this workshop by a package of subroutines called PGPLOT. A manual describing this set of subroutines and how to use them can be borrowed from the Laboratory Preparation Room or you can use the links on the World Wide Web page:
http://www.star.le.ac.uk/~rw/compshop/
REAL(SING_PREC), ALLOCATABLE, DIMENSION(:) :: X,Y,VX,VY,T,PE,KE
This is a declaration of allocatable arrays which must appear before the executable statements. You will notice that the dimension of the arrays is not specified (:). In our program we do not know the number of time steps required until it is typed in by the user therefore we don't know beforehand how large the arrays need to be. After the line which reads in NSTEP we can actually allocate the array space using:
ALLOCATE(X(NSTEP),Y(NSTEP),VX(NSTEP),VY(NSTEP),T(NSTEP))
ALLOCATE(PE(NSTEP),KE(NSTEP))
The loop which performs the calculation should then be modified as follows:
TIME=0.0
X(1)=P(1)
Y(1)=P(2)
VX(1)=P(3)
VY(1)=P(4)
T(1)=TIME
PE(1)=PEND_PE(P)
KE(1)=PEND_KE(P)
DO J=2,NSTEP
CALL RUNKUT(DT,PEND_EQM,TIME,NV,P)
X(J)=P(1)
Y(J)=P(2)
VX(J)=P(3)
VY(J)=P(4)
T(J)=TIME
PE(J)=PEND_PE(P)
KE(J)=PEND_KE(P)
ENDDO
The PGPLOT calls needed to plot the results are:
! Open plotting device
CALL PGBEGIN(0,"?",1,1)
! Define environment (XMIN and XMAX must be declared SING_PREC)
XMIN=MINVAL(X)
XMAX=MAXVAL(X)
XRANG=XMAX-XMIN
XMIN=XMIN-XRANG*0.05
XMAX=XMAX+XRANG*0.05
CALL PGENV(T(1),T(NSTEP),XMIN,XMAX,0,1)
! Draw axes and border with labels
CALL PGLABEL("time","X","Pendulum displacement vs. time")
! Draw curve
CALL PGLINE(NSTEP,T,X)
! Clear up and end plotting
CALL PGEND
The library of subroutines must be made available to your program when you use the compiler/linker and the commands you will need to produce a plot are as follows.
$ . /disk/s/zrw/under/pgplot_f90 $ f90 <name>.f90 `pgplot_f90_link` -o <name> $ <name>
The first line (starting with dot space) is only required once during a session. The `pgplot_f90_link` (use the single quotes at the top left of the keyboard) includes the necessary object libraries into the f90 command. When you run the program the PGBEGIN routine will prompt for a device name. The device type you should use on the X-terminals is /XS (or /GF if you are using a serial line terminal). Hardcopy can be obtained using /PS (landscape) and /VPS (portrait). These create a Postscript text file called pgplot.ps which may be sent to the print queue using the following command:
$ lp -dA4 < pgplot.ps
You can get the hardcopy from the Print Station in the terminal room between the 1st and 2nd year Laboratories. You will need a Print Copycard which you can get from the demonstrator.
A few more words about arrays are in order. In the functions PEND_PE and PEND_KE the dimension of the array P(4) was declared in both the calling program and the function and there was no mention of the size in the argument list. In this context, where the number of elements was fixed by the problem (the dimension of the phase space of the pendulum), this was reasonable. However in general it is a good idea to pass the size of array dimensions as INTEGER arguments to the FUNCTION or SUBROUTINE. The start of the RUNKUT routine illustrates this.
SUBROUTINE RUNKUT(H,F,X,N,Y)
! Fourth order Runge-Kutta
IMPLICIT NONE
INTEGER, PARAMETER :: DOUB_PREC=SELECTED_REAL_KIND(P=10,R=30)
INTEGER :: N ! number of y's
REAL(DOUB_PREC) :: H,X,Y(N) ! x step size, x, y's
EXTERNAL F ! subroutine F(X,Y,YPRIME)
! F returns derivatives wrt x in yprime
! Author Dick Willingale 1996-Jul-25
! Local automatic arrays for workspace
REAL(DOUB_PREC) :: F1(N),F2(N),F3(N),F4(N)
REAL(DOUB_PREC) :: Y1(N),Y2(N),Y3(N),Y4(N)
REAL(DOUB_PREC) :: K1(N),K2(N),K3(N),K4(N)
REAL(DOUB_PREC) :: X1,X2,X3,X4
INTEGER :: J ! local index
...
The N is passed as an argument along with Y and the declaration is of the form Y(N). These arguments are called dummy arguments. The actual location of N and Y is in the calling routine not within RUNKUT. The advantage of this is that the routine can be used for any dimension of y (number of variables in the equations), not just the problem in hand for which N=4 and of course N is available for use within the routine as well.
This routine also demonstrates another useful way of declaring arrays in SUBROUTINEs and FUNCTIONs. In the last few declaration lines a number of local arrays are declared using the argument N (e.g. F1(N)). These are called automatic arrays and can be declared in this manner to create workspace for the calculations inside a SUBROUTINE or FUNCTION.
There is another way of controlling the flow of the execution of statements in F90 called the CASE construct. It is designed specifically for use by character strings or integers and works for logicals. It has the form:
SELECT CASE(expression_select)
CASE(expression1)
F90 statements
CASE(expression2)
F90 statements
CASE DEFAULT
F90 statements
END SELECT
The expression select is a character string (or integer) which will determine which of the following cases is executed. If none of the expression1 etc. is satisified then the DEFAULT will be used. The expressions use must not be REAL. To control flow with REAL values they must be converted to INTEGER or LOGICAL (using some expression) and used in CASE or IF constructs.
If you play about with the starting point of the pendulum bob (initial position and velocity) you will soon discover that the motion of the pendulum can be very complicated with no apparent structure or regular behaviour. It is very difficult to make sense of what is going on by plotting the path as described above.
A more useful view of the motion is provided by what is known as the Poincaré Section, particularly if the air drag (energy loss) is zero (d=0). If you define a plane in phase space, for example the y=0 plane, then, as the pendulum swings about, the path in phase space will cross the plane on a regular basis. Every time the path crosses the plane it defines a point (which can be plotted). Actually we should restrict the points to those for which the path is travelling from y<0 to y> 0 (or the other way around). The section of the path between such points then defines what is called an orbit in phase space. If the motion is allowed to continue then the plotted points gradually build up a picture of the Poincaré Section of the motion.
If the motion is periodic then the path will eventually return to a point already plotted and the sequence will start again. Indeed, if the motion is particularly simple then the picture may consist of just one point! If the motion is chaotic then the bob never returns to the same position in phase space and the Poincaré Section picture will consist of an infinite number of points covering some area of the plane. Somewhere in between the simple regular picture and chaos there are some rather complicated patterns as you will see.
The next task is to modify your program to plot a Poincaré Section. The changes are remarkably simple. You must now control an outer loop by the number of section points required rather than the number of steps. An inner loop must be indefinite since you no longer know how many integration steps are required. If the number of points requested is NP then the loops will look something like:
TIME=0.0
DO J=1,NP
DO
POLD=P
CALL RUNKUT(TINC,PEND_EQM,TIME,NV,P)
IF(P(4)>0.0.AND.POLD(2)<0.0.AND.P(2)>=0.0) THEN
EXIT
ENDIF
ENDDO
X(J)=POLD(1)
Y(J)=POLD(2)
VX(J)=POLD(3)
VY(J)=POLD(4)
PE(J)=PEND_PE(POLD)
KE(J)=PEND_KE(POLD)
T(J)=TIME
ENDDO
An array POLD(4) is used to save the position before RUNKUT is called. The inner loop keeps running until the path crosses the y=0 plane going in the positive direction. Since the motion of the pendulum is bounded this is going to happen sooner or later. The outer loop saves NP points of the section.
You can plot out the points using:
CALL PGPOINT(NP,X,VX,ICHAR)
where ICHAR=1 (type INTEGER) gives a point. To plot other symbols use higher values of ICHAR. For example ICHAR=2 gives + and ICHAR=3 gives *.
The first and second year Fortran Workshops have introduced most of the important features of F90. However there are a few topics which have been omitted. F90 can handle pointers to data rather than using conventional variables. These are most useful when the data base is dynamic, for example in real-time control and analysis. The power and scope of I/O has been ignored. In particular we have not mentioned unformatted input/output and have only touched apon FORMAT statements for contolling formatted I/O. There are a very large number of intrinsic procedures (FUNCTIONs like SIN, MAXVAL etc.) which provide many facilities. Before embarking on writing a new program you should check that F90 doesn't already provide a way of doing exactly what you want.
26 letters of the alphabet A-Z uppercase only (Fortran compilers ignore case) 10 decimal digits 0-9 4 arithmetic symbols + - * / (note ** for raise to the power) 18 other symbols ( ) . = , ' $ : space ! " % & ; < > ?
; semi-colon is used to separate statements on a line ! exclamation indicates the start of comment (rest of line ignored) & ampersand at end of line indicates continuation on next line
Integer constants 1 411 -34 +18 Real constants 1.35 44.33 -0.0123 1.3E-3 Character strings "This is a character string"
variable names up to 31 characters alphabetic, underscore or numeric first character must be alphabetic without IMPLICIT NONE or some other IMPLICIT declaration then: if first character I,J,K,L,M,N then type INTEGER otherwise type REAL
All declarations within a code module must appear before any executable statements.
IMPLICIT NONE ! no name convension for type
INTEGER :: I,IFLAG ! declare integer variables
INTEGER, PARAMETER :: INUM=100 ! declare integer parameter (constant)
INTEGER, PARAMETER :: D_P=SELECTED_REAL_KIND(P=10,R=30)
! define precision parameter
REAL :: RR, ZZ ! declare real variables
REAL(D_P) :: DD,FF ! real variables with precision
REAL, DIMENSION(3) :: POS ! array dimension
REAL :: TIM(3) ! array dimension declaration
LOGICAL :: FLAG ! logical variable (true or false)
CHARACTER(LEN=40) :: BLURB ! character variable
CHARACTER(LEN=10), DIMENSION(3) :: TARR
! character array of 3 elements
! each element 10 characters
REAL, ALLOCATABLE, DIMENSION(:) :: X,Y
! dynamically allocatable arrays
REAL :: A=1.45 ! declare a real with initial value
REAL, DIMENSION(2) :: B=(/ 1.23,3.45 /)
! initial array element values
REAL, DIMENSION(50) :: C=(/ (0,I=1,50) /)
REAL, DIMENSION(40) :: D=0.0 ! implied/implicit loops for initial values
The actual allocation of memory for dynamically allocatable arrays must be done within the executable code after the declarations:
... NSTEP=1000 ... ALLOCATE(X(NSTEP),Y(NSTEP)) ! allocate array space as executable ...
I=?? ! get integer subscript from somewhere X(I)=5.0 ! reference array element using integer Y(I)=X(I)+10.0
READ *, R1,R2 ! use the default input stream PRINT *, "The sum is ",SUM ! use the default output stream
OPEN(UNIT=1,FILE="flight.dat",STATUS="OLD") ! open an existing file READ(UNIT=1,FMT=*) T,H ! read from a unit number ... OPEN(UNIT=2,FILE="output.dat",STATUS="NEW") ! open a new file WRITE(UNIT=2,FMT=*) T,H ! write to a unit number ... CLOSE(UNIT=1) ! close a unit number
PRINT 10, II,RR ! format specified by a labelled statement
10 FORMAT(5X,I3,1X,F10.3)
WRITE(2,"(5X,I3,1X,F10.3)") II,RR ! format specified in-line
The following field descriptors are available:
nX n spaces
Fn.m REAL number, fixed point, m digits after point
En.m REAL number with exponent, m digits after the point
In INTEGER
An CHARACTER string
Arithmetic operators:
( ) parentheses to dictate the order of operations ** exponentiation ( raising to a power) * or / multiplication or division + or - addition or subtraction
INTEGER :: J REAL : X,Y,A,B ... X=A+B ! set variable X equal to A+B J=J-1 ! decrement variable J by 1 Y=X**2 ! set variable Y to X squared
INTEGER :: I,L,II REAL : A,R,R1,R2 ... I=R+35.3 ! The REAL result is truncated before assignment to INTEGER A=I/L ! The result is truncated to an INTEGER and assigned to REAL R2=R1+II ! The result of mixed arithmetic is REAL
The following relational operators are available:
== or .EQ. equals /= or .NE. not equals > or .GT. greater than < or .LT. less than >= or .GE. greater or equal <= or .LE. less or equal
The following logical operators are available:
.NOT. invert logical value .AND. logical and .OR. logical or .EQV. logical equivalence .NEQV. logical non-equivalence
Decision structures:
IF(A<3.0) THEN ! start of an IF structure ... ELSE IF(A>10.0) THEN ! continue IF structure ... ELSE ! final clause of IF structure ... END IF ! end of IF structure
SELECT CASE(expression_select) ! start of a CASE structure
CASE(expression1) ! expression CHARACTER or INTEGER
! do something
CASE(expression2)
! do something else
CASE DEFAULT
! some default
END SELECT ! end of CASE structure
Definite loops:
DO J=1,10 ! start a definite loop with counter J
! within loop J will be 1,2,....9,10
END DO
DO K=1,13,2 ! start definite loop with increment 2
! within loop K will be 1,3,5....9,11,13
END DO
Indefinite loops:
NMAX=100 ! set upper limit
NMIN=-100 ! set lower limit
DO ! start indefinite loop
J=?? ! get J from somewhere, somehow!
IF((J>NMAX).OR.(J<NMIN)) THEN ! test for exit
EXIT ! jump to after END DO
ELSE
! keep going
ENDIF
END DO
...
Implied loops:
REAL , DIMENSION(10) :: A ... READ(UNIT=1,FMT=*) (A(K),K=1,10) ! implied loop on input
implicit loops for arrays:
REAL, DIMENSION(100) :: A,B,C ... A=1.4 ! set all elements of array A to a scalar B=A+4.0 ! add 4.0 to all elements of A C=A*B ! multiply elements of A by elements of B D=SIN(C) ! take the SIN of each element of C
ABS absolute value SQRT square root SIN sine (argument radians) COS cosine (argument radians) TAN tangent (argument radians) ASIN arcsine (returns radians) ACOS arccosine (returns radians) ATAN arctangent (returns radians) ATAN2(X,Y) arctangent Y/X (returns radians) ALOG10 logarithm base 10 EXP raise e to power INT truncate REAL to INTEGER NINT fix REAL to nearest INTEGER REAL float INTEGER to a REAL MOD remainder of A1/A2
A single PROGRAM module provides the entry point from the operating system:
PROGRAM progname ... END PROGRAM progname
SUBROUTINEs can be called from within PROGRAM or other SUBROUTINEs or FUNCTIONs:
SUBROUTINE subname(dummy argument list) ... END SUBROUTINE subname
FUNCTIONs can be used within PROGRAM, SUBROUTINEs, or other FUNCTIONs:
FUNCTION funname(dummy argument list) ... END FUNCTION subname
A MODULE is used for the declaration of variables and parameters that are to be shared (global) between PROGRAM, SUBROUTINEs and FUNCTIONs:
MODULE modname
declarations...
END MODULE modname
SUBROUTINE whatever
USE modname ! use declarations in a MODULE
...
Calling subroutines and using functions:
SUBROUTINE DOIT(A,B)
IMPLICIT NONE
REAL :: A,B ! declare dummy arguments
internal declarations
executable instructions
END SUBROUTINE DOIT
FUNCTION ZAP(C)
IMPLICIT NONE
REAL :: ZAP,C !declare function and dummy arguments
internal declarations
executable instructions
ZAP=expression !assign value to function before return
END FUNCTION ZAP
PROGRAM TRY
REAL :: X,Y,Z
X=4.0
CALL DOIT(X,Y)
Z=ZAP(X)
PRINT *, X,Y,Z
END PROGRAM TRY
A safe way of passing arrays as arguments:
SUBROUTINE THINGY(N,A,B)
IMPLICIT NONE
INTEGER :: N
REAL, DIMENSION(N) :: A
REAL :: B
...
END SUBROUTINE THINGY
PROGRAM TRYIT
INTEGER, PARAMETER :: NEL=10
REAL :: X(NEL),Y(NEL)
...
CALL THINGY(NEL,X,Y)
...
END PROGRAM TRYIT