-------------------------- WHY LEARN FORTRAN? ------------------------ Most IGPP faculty program in FORTRAN. For many years, FORTRAN was the language of choice for scientific programming. Recently, C has emerged as a competitor and seems to now be much more widely taught to students, perhaps because it is favored by computer science departments (who worry more about writing compilers than how to handle complex numbers). Why learn FORTRAN if you already know C? Several reasons come to mind: (1) You will communicate and exchange software more readily with most IGPP faculty, who are generally proficient in FORTRAN but rather challenged in C. (2) There is a huge library of existing FORTRAN subroutines to perform various algorithms, both at IGPP and in the wider scientific community. (3) Complex numbers and double precision are included. (4) It's fun! Finally, some cranky advice: It is important to become familiar with at least one major programming language. FORTRAN and C are such languages; MATLAB and other application programs are not. In the long run, you will handicap your ability to do science if you do not take the time to gain experience with a real programming language. In the worst case, you will be one of those annoying people that is forever asking others to write or modify programs for you. ------------------------ FORTRAN HISTORY ----------------------------- The name Fortran is derived from IBS Mathematical FORmula TRANslation System. Originally it was spelled in all caps as FORTRAN, but the more modern usage is to only capitalize the first letter. Fortran is one of the oldest computer programming languages and was begun in 1954 by IBM programmers led by John Backus (no relation to George). It is updated and "improved" by some committee of computer sciences every ten years or so. Important versions include: Fortran IV (released in 1972) Fortran 77 (released in 1980) Fortran 90 (released in 1991) Fortran 95 (released in ?, not installed on our Suns) Most versions are designed to be fully backward compatible with previous versions. However, Fortran90 departed from the column sensitive format of older Fortran, resulting in a more modern approach that necessitated some slight incompatibilities with older Fortran. I will try to teach this class entirely in Fortran90, but will describe the differences with Fortran77 when they are important because you are likely to need to work with existing Fortran77 code at some point. ------------------------ FORTRAN on the IGPP Macs -------------------- All of the Barnyard Macs should have F77 installed because it can be downloaded for free. Unfortunately there is not yet a good free F90 compiler for the Macs. We will be using the IBM F90 compiler, which costs $400 per machine. We are looking into getting a site license but have not been able to get a significant discount yet. Accordingly, we have installed the F90 compiler only on grunt and you will have to login to grunt to compile your programs. Remember that you can do this remotely by entering: ssh -X grunt.ucsd.edu Please let me or Graham know immediately if you have any problems getting the F90 compiler to work! Many of the comments in these notes still refer to the Suns (and my old desktop machine, rock) but all of the programs should work equally well on the Macs. As we go through the notes we will learn if there are any Mac compatibility issues. You can edit your Fortran source code with any plain ascii text editor (don't try to use Word!) but on the Macs I recommend that you use the XCode program which has lots of nice features for programming. To use Xcode conveniently, I have put the following line into my .cshrc file: alias xcode 'open /Developer/Applications/Xcode.app/' With this alias, I can just enter xcode progname.f90 *****Aaakkkk! So far I can't get Xcode to work with the .f90 suffix, just with the .f suffix. Perhaps somebody can do some research on this to find out if there is some setting to change--we need to be able to use it with .f90 but for now we will have to use a regular text editor. ---------------------- TEXT and MANUALS --------------------------- This class will be a tutorial on Fortran90 and will not be comprehensive. Thus, I recommend that you buy a textbook that will be serve as a more complete reference. I have been using Guy's copy of: Programmer's Guide to Fortran 90, Brainerd, Goldbert and Adams, Springer, 1996. which looks pretty good. It appears to be going out of print, but Amazon.com lists a few new copies for around $60. A newer version of this book is (I think): Fortran 95 Handbook, Adams (editor) , Brainerd, Martin and Weqener, MIT Press, 1997. It has 750 pages and lists for $65. Amazon also lists: Fortran 90/95 Explained, Metcalf and Ker Reid, Oxford Univ. Press, 1999 which has some favorable reviews. It's only $40 so I will probably buy it. The UCSD bookstore has only one Fortran book: FORTRAN 90 for Engineers and Scientists, Nyhoff and Leestma, Prentice Hall, 1996. It is long (1071 pages) and expensive ($76) and I don't really recommend it. You may want to check around to see which book you like best. Online manuals for the Sun Fortran90 compiler may be obtained at http://docs.sun.com (search on Fortran 90 to find the Handbook and the Users Guide). ---------------- HOW TO COMPILE FORTRAN90 PROGRAMS ------------------- FORTRAN90 programs are written as ascii files that end in ".f90". This is called the source code. Programs in older versions of Fortran end in simply ".f" Because of the slight incompatibilities between the versions, be sure to use the appropriate suffix so that both you and the compiler know which version to use. It is in fact possible to set things up so that your F90 programs end in .f rather than .f90. I do not recommend this because there is so much existing code around in Fortran77 that confusion is likely to result if you do not explicitly identify your code as F90. Before the program can be run, it must be compiled using the FORTRAN90 compiler on your computer, creating the executable file that you use to actually run the program. For example, suppose you have a program called printmess.f90 with contents: ! simple Fortran90 test program (printmess.f90) program printmess print *,"test message" end program printmess To compile this program from my computer (rock), I enter "f90 printmess.f90 -o printmess" and get the following response: rock% f90 printmess.f90 -o printmess rock% No news is good news! If there was a syntax or other error detected during the compilation, we would get an error message at this point. This creates an executable version of the program called "printmess" (you should ALWAYS use name of the program without the ".f90" for the executable) which can then be run simply by entering printmess: rock% printmess test message rock% The "f90" command invokes the Sun FORTRAN90 compiler. "f77" is the corresponding command for FORTRAN77. If you are like me, you will quickly tire of typing in the line "f90 printmess.f90 -o printmess" Thus, I recommend that you create a Makefile (yes, the first character MUST be capitalized) in the same directory as the program. Include in the Makefile the following lines: %: %.f90 f90 $< -o $* ***The tricky part of this is that the space before the f90 must be entered as a tab, not as a series of spaces.*** Once you have set up this Makefile, then you can just enter: make printmess in order to compile the program. Makefiles are very useful to keep track of compiler options and to bind with subroutines. (more about this later). The difference between source code and excutables is fundamental to languages like FORTRAN and C. It is why they generally run much faster than uncompiled languages like "BASIC" or MATLAB scripts. ---------------------- THE FIRST PROGRAM EXPLAINED ------------------- OK, now let's examine our simple F90 test program again: ! simple Fortran90 test program (printmess.f90) program printmess print *,"test message" end program printmess The first line is a comment line. Anything following an exclamation mark (!) is a comment. The end of the line terminates the comment. There is no need (as in many languages) to terminate the comment with another flag. We can also use "!" to add an inline comment following a statement, e.g., print *,"test message" !print message on screen would be a valid line. The next line is blank. You can put in as many blank lines as you want and they will be ignored by the compiler. Blank lines provide a good way to improve the readability of longer programs by breaking them up into coherent blocks of code. The next line is used to name the program. The name printmess is not used by the program at all. This line is optional but is considered good programming style. Good style also suggests that lines in the body of the program are indented, in this case by three spaces: print *,"test message" print * will output to the screen. The desired output is enclosed in quotes (apostrophes would also work). Finally, all F90 programs must terminate with an "end" statement. The "program printmess" is optional, but is considered good programming practice. These style matters don't matter much for a short program like this, which would work just as well if were written as: print *,"test message" end but are helpful for longer programs (hundreds to thousands of lines of code) where the label following the "end" would remind a reader of which program is actually ending. ASSIGN F1 Write a F90 program to print your favorite pithy phrase. --------------------- HOW TO MULTIPLY TWO INTEGERS ------------------ Here is a simple F90 program to multiply 2 and 3: program multint integer :: a, b, c !declare variables a = 2 b = 3 c = a * b print *,"Product = ",c end program multint The program uses three variables, the letters a, b and c. Variable names can be from 1 to 31 characters long. The first character must be a letter. The remaining characters can be any combination of letters, numbers, and underscores (_). Variables in Fortran can be of many types, including real (floating point), integer, complex, double precision, character and logical. In our case, we want them to be integers so we define them using a "type statement" integer :: a, b, c !declare variables (older Fortran programs do not include the :: in these statements; this convention still works under F90 but is discouraged (why?!) ) For the Sun F90 compiler, these are 4-byte integers that can range from -2,147,483,648 to 2,147,483,647. Alternatively, they could have been defined as real numbers: real :: a, b, c in which case they would be floating point (real) variables. For the Sun compiler, real numbers range from 1.175494e-38 to 3.402823e+38) The lines that follow are pretty self explanatory: a = 2 b = 3 c = a * b print *,"Product = ",c Note that * is used to indicate multiplication as is true in almost all programming languages. Addition is +, subtraction is -, and division is /. In Fortran, a to the power of b is written as a**b where a and b can both be real, another Fortran advantage over standard C. ---------------- AN IMPORTANT HISTORICAL DIGRESSION ------------- It is not required that variables be declared in Fortran. If they are not declared, then the compiler assigns them as integer or real, depending upon their first letter. Variable names beginning with i through n (INteger, get it?) are assumed to be integers; all others are assumed to be real. Many, if not most, older Fortran programs adopt this convention. Often they only declare variables when the rule is broken, for example when it is desired that the variable "year" be an integer. A prominent example of this type of Fortran programming may be found in the first edition (1986) of Numerical Recipes (but by the time of the second edition in 1992, they had been shamed into declaring all variables). Such undeclared variables are said to be declared implicitly based upon their first letter. An "implicit" statement can be used to modify the default rules, for example implicit real (a-z) will make all undeclared variables real. However, modifications such as this will only lead to more confusing code. Modern programming practice is to explicitly declare ALL variables and to have the compiler warn us when variables are not declared. This can be done by including the staement implicit none at the beginning of every Fortran program. I will admit that, as a long time Fortran programmer, I do not always follow this practice. I must concede, however, that it is a good idea and is likely to save more time in the long run (by eliminating program bugs that are often caused by undeclared variables) than is lost while writing the program. Example 1: Suppose you should have the following line in your program: x2 = a1 + a2*sin(theta)*scale1/scale2 but you accidentally type: x2 = a1 + a2*sin(theta)*scalc1/scale2 If you don't follow the practice of declaring all of your variables, then the error will not be detected during compilation. The program will run and assign zero to the otherwise unused variable "scalc1" and you will get wrong answers. On the other hand, if you use "implicit none" and declare your variables, the compiler will flag "scalc1" as an undeclared variable and you can fix it before it causes any more trouble. Example 2: Suppose you have the following lines in your code: kmdeg = 111.19 dist = (delta2 - delta1)*kmdeg but you have not declared kmdeg explicitly. Because the letter "k" is between i and n, the program will declare kmdeg as an integer. The value 111.19 will be truncated to 111 and you will get values for dist that are slightly, but not obviously, wrong (the worst kind of program bug to have). To save you from these embarassments and to encourage good programming habits, we will declare all variables for the programs in this class. ASSIGN F2 Copy the program multint but leave out the integer :: a, b, c statement. What happens when you run the program? Why? ASSIGN F3 Cut and paste the following defective program onto your computer: program longjump beamon_long = 8.90 !distance in meters powell_long = 8.95 dif_inch = (powell_long - beamon_1ong) * 39.37 print *,"Powell jumped ",dif_inch," inches more than Beamon" end program longjump Compile and run the program. Then explain why it gives the wrong answer. ------------------ ALTERNATIVE CODING OPTIONS ------------------ There are always lots of ways to write the same program. Here is another way to write the multint.f90 code: program multint2 implicit none integer :: a=2, b=3, c !declare variables c = a * b; print *,"Product = ",c end program multint2 First, notice that variables can be assigned values when they are declared (OK in F90, don't try this in F77). Second, notice that more than one command can be included on a line if they commands are separated by a semicolon (again, only OK in F90). Both of these changes make the code more similar to C. The variable assigning option is a reasonable convenience, but the multiple command option is certainly misused in this case because it makes the code much harder to read. Unless you have a really good reason to put more than one command on a line (saving space is NOT a good reason!), I suggest that you never use semicolons in this way. --------------- MAKING A TRIG TABLE USING A DO STATEMENT -------------- Any reasonable programming language must provide a way to loop over a series of values for a variable. In C, this is most naturally implemented with the "for" statement. In FORTRAN this is done with the "do" loop. Here is an example program that generates a table of trig functions: program trigtable implicit none real theta, stheta, ctheta, ttheta, degrad degrad = 180./3.1415927 do theta = 0.0, 89.1, 1.0 ctheta = cos(theta/degrad) stheta = sin(theta/degrad) ttheta = tan(theta/degrad) print "(f5.1,1x,f6.4,1x,f6.4,1x,f7.4)", theta, ctheta, stheta, ttheta enddo end program trigtable ****IMPORTANT NOTE**** This example uses a real do variable. This is permitted in F77 but is NOT standard F90. This example will compile with the Sun and IBM F90 compilers but will give an error with the Mac g95 compiler. So I need to change this example to not use a real do variable. Fortran (like C and Matlab) uses radians (not degrees) as the arguments of the trig functions. Thus, following the definitions of the variables as real, we assign degrad to 180/pi so that we can easily make this conversion. do theta = 0.0, 89.1, 1.0 This begins the "do" loop which must eventually be closed with the "enddo" statement. For clarity, we indent the inside of the loop. This is a loop over values of theta from a starting value of 0.0, incremented by 1.0 each time, until theta is greater than 89.1. Thus theta will assume the values (0.0, 1.0, 2.0, ...., 88.0, 89.0) inside the loop. Note that we use 89.1 rather than 89.0 as the ending value to avoid the possibility that roundoff error might cause the desired ending value (computed by successively adding 1.0 to theta) to slightly exceed 89 and thus be excluded from the loop. These roundoff peculiarities prompt some books to advise against using real arguments in do loops. However, I find them too much of a convenience to give up. Inside the do loop we compute the cosine, sine and tangent of theta, after converting from degrees to radians by dividing by degrad. To make the output look nice, we do not use print *, theta, ctheta, stheta, ttheta which would space the numbers irregularly among the columns. Instead, we explicitly specify the output format using a format specification: print "(f5.1,1x,f6.4,1x,f6.4,1x,f7.4)", theta, ctheta, stheta, ttheta where f5.1 specifies that theta will be output as a real number into 5 total spaces, with 1 digit to the right of the decimal place. Similarly, f6.4 specifies that ctheta is output into spaces with 4 digits to the right of the decimal place. The numbers will be right justified, with leading blanks used as necessary. The 1x specifies that one blank character will be output between each of the numbers. An alternative way to write this: print "( f5.1, f7.4, f7.4, f8.4)", & theta, ctheta, stheta, ttheta Here we have used the continuation character & to split the statement into two lines. Blanks are ignored so we can space things out neatly to line up the variables with their formats. Notice that we have removed the need for the 1x between formats by adding an additional column to the appropriate formats (i.e., writing f7.4 rather than f6.4). Because the numbers are right justified, this will add an additional space to the left of each number. Aligning things this neatly is probably more trouble than it's worth, but it certainly makes the code easier to understand. Notice that in this case the f7.4 appears twice in a row. Often programmers will write this more compactly as: print "(f5.1, 2f7.4, f8.4)", theta, ctheta, stheta, ttheta since Fortran allows this syntax. Older Fortran programs usually put the format specifier into a separate numbered line, i.e. print 117, theta, ctheta, stheta, ttheta 117 format (f5.1, 2f7.4, f8.4) This convention is still allowed in F90 but should not be used unless you want to use the format statement more than once, e.g., print 117, theta1, ctheta1, stheta1, ttheta1 print 117, theta2, ctheta2, stheta2, ttheta2 117 format (f5.1, 2f7.4, f8.4) 117 is termed a line label and must consist of digits. It is best not to refer to it as a line number (the old usage) because it normally has nothing to do with the line numbers and the labels need not be sequencial (more about this later). Note that the format line need not immediately follow the print line(s); it can come before. Many older programs put all of the format statements at the end of the code. An altenative way in F90 to use the same format specification more than once is to define it as a character variable, i.e., character (len = 30) :: fmt fmt = "(f5.1, 2f7.4, f8.4)" print fmt, theta1, ctheta1, stheta1, ttheta1 print fmt, theta2, ctheta2, stheta2, ttheta2 but we are getting ahead of ourselves because we have not yet shown how to use character variables. ------------------ FORTRAN MATHEMATICAL FUNCTIONS--------------------- We used the Fortran sine, cosine and tangent function in the trigtable program. Here is a complete list of math functions: acos(x) arccosine asin(x) arcsine atan(x) arctangent atan2(y,x) arctangent of y/x in correct quadrant (***very useful!) cos(x) cosine cosh(x) hyperbolic cosine exp(x) exponential log(x) natural logarithm log10(x) base 10 log sin(x) sine sinh(x) hyperbolic sine sqrt(x) square root tan(x) tangent tanh(x) hyperbolic tangent ------------------ POSSIBLE INTEGER/REAL PROBLEMS -------------------- As an aside, note that the trigtable program uses: degrad = 180./3.1415927 rather than simply degrad = 180/3.1415927 The reason is to make completely sure that the program will compute a real quotient and not an integer quotient. In fact, this caution is not needed in this case, as the following program demonstrates: program testfrac implicit none real c c = 2/3 print *,'2/3 = ', c c = 2/3. print *,'2/3. = ', c c = 2./3 print *,'2./3 = ', c c = 2./3. print *,'2./3. = ', c end program testfrac Running the program yields: 2/3 = 0.0E+0 2/3. = 0.6666667 2./3 = 0.6666667 2./3. = 0.6666667 As long as one part of the fraction is real, the program will compute a real quotient. It is only when both numbers are written as integers that the result is truncated. However, I have gotten into the habit of always including the decimal point in real expressions to avoid someday accidentally writing something like: a = sin(phi/degrad) + (2/3) * cos(theta/degrad)**2 which will definitely produce the wrong answer! --------------------- MORE ABOUT FORMATS ----------------------- There are many different formats that can specified. Here are some common examples: i5 = integer, 5 spaces, right justified i5.4 = as above but pad with zeros to 4 spaces (88 becomes 0088) f8.3 = real, 8 spaces, 3 to right of decimal place e12.4 = real output with exponent, 4 places to right of decimal e.g., b-0.2342E+02 where "b" is blank (useful for big or small numbers or when you are not sure what size they will be and want to be sure you have enough room to output them) If a number does not fit into the allocated spaces, it will appear as a series of astericks (*****) a8 = character output, 8 spaces, right justified (if the length of the string is greater than 8, then the leftmost 8 characters will appear) 2x = output two blanks tn = tab to position n tln = tab left n spaces trn = tab right n spaces (tr2 is the same as 2x) / = start a new line ASSIGN F4 Write a F90 program to print a table of x, sinh(x), and cosh(x) (the hyperbolic sine and cosine, these are built-in functions in Fortran) for values of x ranging from 0.0 to 6.0 at increments of 0.5. Use a suitable format to make nicely aligned columns of numbers. ------------------ INPUT USING KEYBOARD -------------- So far all of our example programs have run without prompting the user for any information. To expand our abilities, let's learn how to input data from the keyboard. In most programs, we will want to first prompt the user to input the data, so here is an example of how to input two numbers: print *, "Enter two integers" read *, a, b Pretty easy, isn't it? (Compare this section with the corresponding input section in the C notes and you will see why I think Fortran is more user friendly than C) Here is an example of a complete program that multiplies two numbers: program usermult implicit none integer :: a, b, c print *, "Enter two integers" read *, a, b c = a * b print *, "Product = ", c end program usermult Running this program, we have: rock% usermult Enter two integers 2 5 Product = 10 The program will also accept the input on two lines: rock% usermult Enter two integers 2 5 Product = 10 This often happens to me if I forget that I'm supposed to input more than one number. When the program just sits there, I then realize that I need to input more numbers (or I wonder why it's taking so long to finish!). What happens if we make a mistake and try entering a real number? Let's check: rock% usermult Enter two integers 3.1 15 ****** FORTRAN RUN-TIME SYSTEM ****** Error 1083: unexpected character in integer value Location: the READ statement at line 4 of "usermult.f90" Unit: * File: standard input Input: 3.1 ^ Abort rock% The error message in this case is quite informative and tells us exactly what the problem is. This is better than performing the computation and returning the wrong answer (the default in C for this example unless you are quite careful). ------------------- IF STATEMENTS ------------------- Next, let's modify this program so that it will allow the user to continue entering numbers until he/she wants to stop: program usermult2 implicit none integer :: a, b, c do print *, "Enter two integers (zeros to stop)" read *, a, b if (a == 0 .and. b == 0) exit c = a * b print *, "Product = ", c enddo end program usermult2 Here the do loop has no arguments and the block of code inside the do loop will be repeatedly executed until an "exit" command is executed. As in our previous example, we indent the block inside the do loop to make the code easier to read. The program will allow the user to continuing entering numbers to be multiplied. When the user wishes to stop the program (in a more elegant way than hitting C!), he/she enters zeros for both arguments. The "if" statement checks for this and exits the loop in this case: if (a == 0 .and. b == 0) exit A list of the relational operators in different languages is as follows: FORTRAN 77 90 C MATLAB meaning .eq. == == == equals .ne. /= != ~= does not equal .lt. < < < less than .le. <= <= <= less than or equal to .gt. > > > greater than .ge. >= >= >= greater than or equal to .and. .and. && & and .or. .or. || | or The F77 syntax will still work under F90 and you are likely to see this in many of the older programs, e.g., if (a .eq. 0 .and. b .eq. 0) exit will also work. These operators can be combined to make complex tests, e.g., if ( (a > b .and. c <= 0) .or. d == 0) z = a There is, of course, an order of operations for these things which I can't remember very well. Look it up in a book if you are unsure or, better, just put in enough parenthesis to make it completely clear to anyone reading your code. One nice aspect of Fortran compared to C is that if you make a mistake and type, for example, if (a = 0) exit you will get an error message during compilation. In C this is a valid statement with a completely different meaning than is intended! ------------------ IF, THEN, ELSE CONSTRUCTS ------------------ In the above example, a single statement is executed when the if condition is true. A more versatile form is as follows: if (logical expression) then (block of code) else if (logical expression) then (block of code) else if (logical expression) then (block of code) . . else (block of code) end if The blocks of code can contain many lines if desired. As many else if statements as required can be used. At most, one block of code will be executed (once one of the "if" tests is satisfied, it does not check the others). The final "else" will be executed if none of the preceding if statements is true. The final "else" is optional. Here is a demonstration program that repeatedly prompt the user for a positive real number. If it is negative, ask the user to try again. If it is positive, it computes and displays the square root using the sqrt() function. If the user enters zero, the program stops. program usersqrt implicit none real :: a, b do print *, 'Enter positive real number (0 to stop)' read *, a if (a < 0) then print *,'This number is negative!' cycle else if (a == 0) then exit else b = sqrt(a) print *, 'sqrt = ', b end if end do end program usersqrt Notice the use of the "cycle" command, which directs the program to the next iteration of the do loop. In contrast, the "exit" command exits the do loop entirely. The "cycle" and "exit" commands are new to F90 and permit code to be written that is free of the "go to" statements that would likely have been present in a F77 version of this program ("go to" statements are considered mortal sins by the programming style police). ASSIGN F5 Write a program to repeatedly ask the user for the constants a, b, and c in the quadratic equation a*x**2+b*x+c=0. Using the quadratic formula, have the program identify and compute any real roots. Output the number of real roots and their values. Stop the program if the users enters zeros for all three values. --------------------- GREATEST COMMON FACTOR EXAMPLE ------------------ Here is an example program that uses some of the concepts that we have just learned: ! compute the greatest common factor of two integers program gcf implicit none integer :: a, b, i, imax do print *, 'Enter two integers (zeros to stop)' read *, a, b if (a == 0 .and. b == 0) then exit else do i = 1, min(a,b) if (mod(a,i) == 0 .and. mod(b,i) == 0) imax=i end do print *, "Greatest common factor = ", imax end if enddo end program gcf This is not a particularly efficient algorithm, but it runs plenty fast for small numbers. There are some new things here: min(a,b) computes the minimum of a and b using the intrinsic Fortran "min" function. Naturally there is also a "max" function. "min" and "max" can have more than 2 arguments if desired. mod(a,i) computes the remainder of a divided by i (this is called the modulus). If mod(a,i) is zero, then a is evenly divisible by i. If non-zero, then mod(a,i) has the sign of a. Here are some more useful Fortran functions: abs(a) absolute value sign(a,b) abs(a) with sign of b real(a) conversion to real (F77 float(a) still works) int(a) conversion to integer nint(a) nearest integer ceiling(a) least integer greater than or equal to number floor(a) greatest integer less than or equal to number ASSIGN F6 Modify gfc.f90 to compute the least common multiple of two integers. ------------------------ USER DEFINED FUNCTIONS ---------------------- As the length and complexity of a computer program grows, it is a good strategy to break the problem down into smaller pieces by defining functions or subroutines to perform smaller tasks. This provides several advantages: (1) You can test these pieces individually to see if they work before trying to get the complete program to work. (2) Your code is more modular and easier to understand. (3) It is easier to use parts of the program in a different program. To illustrate how to define your own function, here again is the greatest common factor program: program gcf2 implicit none integer :: a, b, getgcf, gcf do print *, 'Enter two integers (zeros to stop)' read *, a, b if (a == 0 .and. b == 0) then exit else gcf = getgcf(a,b) print *, "Greatest common factor = ", gcf end if enddo end program gcf2 function getgcf(x, y) implicit none integer :: getgcf, x, y, i, z do i = 1, min(x,y) if (mod(x,i) == 0 .and. mod(y,i) == 0) getgcf = i end do end function getgcf We now perform the gcf calculation in the function getgcf. The variables a and b in the main program are the function arguments. They are passed to the function in the statement: gcf = getgcf(a,b) The function getgcf definition begins with: function getgcf(x, y) The variables x and y in the function will assume the values passed to the function by the main program. These arguments must match in number and type (real vs. integer, etc.) with the variables in the main program. Note, however, that they do not need to have the same names. The name of the function is by default the value that will be passed back to the main program. In some cases involving recursive functions (a more complicated topic that we may cover later), it may desirable to have the value returned to the main program be specified by a different variable name. To do this, we can use the optional "result" specifier in the subroutine name: function getgcf(x, y) result(z) implicit none integer :: x, y, i, z do i = 1, min(x,y) if (mod(x,i) == 0 .and. mod(y,i) == 0) z=i end do end function getgcf The result(z) indicates that the result will be passed back to the calling program as the variable z. Thus, gcf in the calling program will assume the value of z in the function. Notice that in this case, we do not need to declare getgcf within the function. ASSIGN F7 Modify your program from F6 to compute the least-common multiple as a user-defined function. ---------------------------- SUBROUTINES ----------------------------- Functions are limited in their usefulness because they are designed to pass only one value back to the calling program. A more general construct is the Fortran subroutine, which allows unlimited numbers of values to be passed to and from the calling program. Here is our first geophysically useful example, a subroutine to compute the distance and azimuth between any two points on the Earth's surface: program userdist implicit none real lat1,lon1,lat2,lon2,del,azi do print *,'Enter 1st point lat,lon' read *,lat1,lon1 print *,'Enter 2nd point lat,lon' read *,lat2,lon2 call SPH_AZI(lat1,lon1,lat2,lon2,del,azi) print *,'del, azi = ', del, azi end do end program userdist ! SPH_AZI computes distance and azimuth between two points on sphere ! ! Inputs: flat1 = latitude of first point (degrees) ! flon2 = longitude of first point (degrees) ! flat2 = latitude of second point (degrees) ! flon2 = longitude of second point (degrees) ! Returns: del = angular separation between points (degrees) ! azi = azimuth at 1st point to 2nd point, from N (deg.) ! ! Notes: ! ! (1) applies to geocentric not geographic lat,lon on Earth ! ! (2) This routine is inaccurate for del less than about 0.5 degrees. ! For greater accuracy, use double precision or perform a separate ! calculation for close ranges using Cartesian geometry. ! subroutine SPH_AZI(flat1,flon1,flat2,flon2,del,azi) implicit none real :: flat1,flon1,flat2,flon2,del,azi,pi,raddeg,theta1,theta2, & phi1,phi2,stheta1,stheta2,ctheta1,ctheta2, & sang,cang,ang,caz,saz,az if ( (flat1 == flat2 .and. flon1 == flon2) .or. & (flat1 == 90. .and. flat2 == 90.) .or. & (flat1 == -90. .and. flat2 == -90.) ) then del=0. azi=0. return end if pi=3.141592654 raddeg=pi/180. theta1=(90.-flat1)*raddeg theta2=(90.-flat2)*raddeg phi1=flon1*raddeg phi2=flon2*raddeg stheta1=sin(theta1) stheta2=sin(theta2) ctheta1=cos(theta1) ctheta2=cos(theta2) cang=stheta1*stheta2*cos(phi2-phi1)+ctheta1*ctheta2 ang=acos(cang) del=ang/raddeg sang=sqrt(1.-cang*cang) caz=(ctheta2-ctheta1*cang)/(sang*stheta1) saz=-stheta2*sin(phi1-phi2)/sang az=atan2(saz,caz) azi=az/raddeg if (azi.lt.0.) azi=azi+360. end subroutine SPH_AZI The subroutine is called with the statement: call SPH_AZI(lat1,lon1,lat2,lon2,del,azi) In this case, the lat/lon values are passed to the subroutine while del and azi are passed back to the main program. However, note that if flat1, etc., were changed in the subroutine, then the corresponding variable would also be changed in the main program as well. lat1 in the main program and flat1 in the subroutine point to the same memory location. If one is changed, the other automatically changes as well. Fortran is not case-sensitive so sph_azi and SPH_AZI have the same meaning. I like to put subroutine names in all caps so they are more visible. Note that sph_azi does not have to be declared in the main program. The subroutine is well-documented in this case, explaining exactly what is going into the subroutine and what is going out, as well as some of the limitations of the routine. This may seem like overkill, but documenting your subroutines as completely as possible is likely to save you considerable time later if you ever want to use the routine again. It is good to document the routine well enough that you, or someone else, can use it correctly without having to study the code itself. Clarity is important--I have seen versions of this routine that do not make clear whether the azimuth is measured at the first point to the second point, or vice versa. Listing the limits of the subroutine may help prevent future misuse of the routine, in this example it may prevent the naive user from assuming that the distance returned is accurate when used with geographic latitude and longitude on the Earth. The routine is also designed to be robust with respect to pathological inputs, such as when the two points have the same coordinates. ASSIGNMENT F8 Write a subroutine that computes the volume and circumference of a sphere, given its radius, together with a main program that inputs different values and outputs the results. E-mail me the source code. ---------------LINKING TO SUBROUTINES DURING COMPILATION---------- A powerful aspect of subroutines is that one can link to compiled versions of existing subroutines without having to recompile them. This means that you only have to maintain one version of a subroutine; it need not be listed along with the source code of the main program. For example, the SPH_AZI subroutine is also part of a F77 package of spherical geometry subroutines contained in: /home/shearer/PROG/SUBS/sphere_subs.f Our main program could simply consist of: program userdist2 implicit none real flat1,flon1,flat2,flon2,del,azi do print *,'Enter 1st point lat,lon' read *,flat1,flon1 print *,'Enter 2nd point lat,lon' read *,flat2,flon2 call SPH_AZI(flat1,flon1,flat2,flon2,del,azi) print *,'del, azi = ', del, azi end do end program userdist2 When we compile this program, we need to indicate where the SPH_AZI subroutine can be found. We want to link with what is called an "object file" for the subroutines. Object files end in ".o" and there is a sphere_subs_f90.o file in the same directory (~shearer/PROG/SUBS) as the source file sphere_subs.f. There are slight imcompatibilities between object files generated with F77 versus those generated with F90. To keep track of this, I have started using the naming convention that the F90 object generated from F77 source code has "_f90" just before the ".o" This is why there is both a sphere_subs.o and sphere_subs_f90.o file in my SUBS directory. For my F77 programs, I link to sphere_subs.o while for my F90 programs I link to the sphere_subs_f90.o file. (At first I thought that the object files were compatible between F77 and F90 on the Suns. For many subroutines this turns out to be the case. However, if the routine has a print statement, then they are not compatible and I'm sure there are other examples as well. I have not yet tried this with the IBM XL compiler. Thus it is safer to link only to F90 compiled object files.) To create a F90 object file, enter, for example: f90 -c someprogram.f90 -o someprogram.o You can also create a F90 object file from F77 source code: xlf -c anotherprogram.f -o anotherprogram_f90.o This will work on native F77 code that includes non-F90 compatible features (e.g., comment lines starting with C), because xlf invokes the F77 compiler. Because both xlf and f90 are part of the same IBM compiler package, I assume that they generate compatible code but I have not checked this out completely. It does work for this example but I am curious if it would if there was a print statement in the subroutine (anyone want to try this?). To compile userdist2 and link to sphere_subs.o, we enter: f90 userdist2.f90 /home/shearer/PROG/SUBS/sphere_subs_f90.o -o userdist2 This quickly becomes cumbersome to write so you will find it convenient to set up a Makefile to keep track of all this for you. Here is part of a Makefile that does this for this program: OBJS1= /home/shearer/PROG/SUBS/sphere_subs_f90.o userdist2: Makefile userdist2.f90 $(OBJS1) f90 userdist2.f90 $(OBJS1) -o userdist2 You can list the full path names for any number of subroutine object files in this way. To compile the program, simply enter: make userdist2 ASSIGNMENT F9 Study the SPH_MID subroutine contained in sphere_subs.f (in ~shearer/SUBS). Write a program that uses this subroutine to compute the midpoint between two (lat,lon) points input by the user. Print out the latitude and longitude of this point. Do not attempt to copy the SPH_MID source code, just link to the sphere_subs_f90.o (IBM compiler) or sphere_subs_g95.o (g95 compiler) object file. Make sure that the argument list for SPH_MID in your program matches the subroutine arguments. E-mail me the source code and also tell me where the working executable version of the program is located. Be sure to give me execute permission so that I can try running your program. (Remember: "ls -l" shows the current permissions, "chmod" is how you change the permissions, "man chmod" will explain this.) ----------------------- INTERNAL PROCEDURES ---------------------- The functions and subroutines that we discussed above are called "external" because are located outside of the main program. With external procedures all of the values that are to passed into and out of the procedure must be part of the argument list. All other variables are local to the procedure or to the main program, even if they have the same name as a variable in a different procedure. External procedures are the only kind of procedure allowed in F77 and are what you will find in books like Numerical Recipes or in the subroutine libraries maintained by various scientists at SIO. Their great advantage is their portability--everything you need to know about what they do is contained in their argument list. However, in some cases it may be more convenient to use an "internal" subroutine or function, a method that is new to F90. Internal procedures are listed immediately BEFORE the end statement in the main program. All variable names are shared within internal procedures; thus argument lists are not necessarily required for them to work. Here is an example adapted from the Brainard text that illustrates how internal subroutines work: program sort3 implicit none real :: a1, a2, a3 call read_numbers call sort_numbers call print_numbers contains subroutine read_numbers print *,'Enter three numbers' read *, a1, a2, a3 end subroutine read_numbers subroutine sort_numbers if (a1 > a2) call swap(a1,a2) if (a1 > a3) call swap(a1,a3) if (a2 > a3) call swap(a2,a3) end subroutine sort_numbers subroutine print_numbers print *,'The sorted numbers are: ' print *, a1, a2, a3 end subroutine print_numbers subroutine swap(x,y) real :: x, y, temp temp = x x = y y = temp end subroutine swap end program sort3 Internal procedures are listed following a "contains" statement and before the end statement for the main program. Variables already declared in the main program need not be declared in the internal procedure. Arguments are optional; they are used here in the swap routine to permit it to be used for different pairs from a1, a2 and a3. Note that the procedures are not nested (e.g., one might have tried to put swap internal to sort_numbers); internal procedures may not contain other internal procedures. Note also that variables declared within a subroutine are purely local to the subroutine. For example, the value for temp (in swap) is not available in the main program. Even if temp were declared in the main program, its value will not correspond to the value of temp in the swap subroutine (try it!). The use of internal subroutines is rather pointless in this case because the code the code would probably be clearer without them. However, for longer programs they may well be useful for making the code more structured. If you write a program and notice that you are using the same block of code over and over again, this would a situation where using a subroutine would make sense. The advantage of an internal subroutine is that you don't have to match the argument lists and declare all of the variables. External subroutines can become cumbersome when they involve a large number of variables. Often common blocks are used to avoid long argument lists in this case (more about this later). In many case, internal subroutines may provide a neater way to handle this situation. ------------------------- DOUBLE PRECISION ---------------------- By default, Fortran stores real variables using 4 bytes, providing a precision of about 6 to 7 significant figures over a range from (on the Sun) 1.175494e-38 to 3.402823e+38. If this is insufficient precision, then variables can be defined in double precision. In F77, this was done by declaring them as "double precision" or "real*8" variables. This syntax will still work, although F90 has introduced a new concept--the "kind" specifier for variables. All of these methods are demonstrated in this program: program testdouble implicit none character (len = 30) :: fmt = "(a5,f44.40)" real a4 real*8 a8 double precision aa real (kind=4) :: k4 real (kind=8) :: k8 real (kind=16) :: k16 a4 = 8.9 print fmt, 'a4 = ', a4 a8 = 8.9d0 print fmt, 'a8 = ', a8 aa = 8.9d0 print fmt, 'aa = ', aa k4 = 8.9 print fmt, 'k4 = ', k4 k8 = 8.9_8 print fmt, 'k8 = ', k8 k16 = 8.9_16 print fmt, 'k16= ', k16 end program testdouble which produces the following output on the Suns: a4 = 8.8999996185302734375000000000000000000000 a8 = 8.9000000000000003552713678800500929355621 aa = 8.9000000000000003552713678800500929355621 k4 = 8.8999996185302734375000000000000000000000 k8 = 8.9000000000000003552713678800500929355621 k16= 8.9000000000000000000000000000000003081488 a4 and k4 are regular real*4 variables and approximate 8.9 as 8.999996. Much greater precision is obtained with a8, aa, and k8, which are all real*8 (double precision) variables. k16 is a real*16 variable (quadruple precision). On the Suns, kind=4 is for single precision (4 byte real), kind=8 is for double precision (real*8), and kind=16 is for 16 byte real. Note that the precision of numbers may be specified by appending them with an underscore and the kind value. Thus we write 8.9_8 to indicate that 8.9 is to represented as an 8-byte real number. The use of "kind" is intended to give you greater control over the degree of precision in your programs and ultimately make codes that are less platform dependent in their behavior. I'm not sure if this has been achieved! It appears that F90 also allows complex variables to be declared as double precision, something not allowed in F77. We will check if this is actually true later when we consider complex variables. ***IMPORTANT NOTE: You must define the extra-precision variables using numbers with the appropriate precision. If you write: a8 = 8.9 !***Not correct if a8 is real*8 you will get single precision accuracy. Even the "dble" operator (which works in F77, at least on the Suns) will not work in F90 in assigning variables: a8 = dble(8.9) !***Not correct if a8 is real*8 will assign a8 only at single precision accuracy. The key is to write: a8 = 8.9d0 !correct way to set double precision variables or k8 = 8.9_8 Note that the 0 following the d is the power of 10, thus 1.23d3 is 1230, 0.84d-2 is 0.0084, etc. To get the full 16-byte precision, you must say k16 = 8.9_16 rather than k16=8.9d0 or k16=8.9_8, which will assign k16 only to double precision accuracy. Many of my existing F77 routines use dble( ) to define double precision numbers. These will not work correctly if they are changed to F90, unless all of the dble( ) commands are rewritten. They are fine, however, if they continue to be compiled using F77. I don't know if this is a bug in Sun F90, or if the use of dble( ) in this way is non-standard Fortran. ASSIGNMENT F10 Examine the datetime.f source code in /home/shearer/PROG/SUBS and write a F90 program to compute the number of seconds that have elapsed since noon on your birth date (or the exact time if you know it) until a user specified date and time. Have the program print out this number of seconds. You will want to use the DT_TIMEDIF subroutine. Make sure that all of the variables match and are of the same type (integer, real, and real*8 for timdif, the final argument). You should also be aware that "&" in column 6 of F77 code is how lines are continued, that is: subroutine DT_TIMEDIF(iyr1,imon1,idy1,ihr1,imn1,sc1, & iyr2,imon2,idy2,ihr2,imn2,sc2,timdif) is actually all one line. You should not need to make a copy of these subroutines, just link (as in assignment F9) with the /home/shearer/PROG/SUBS/datetime_f90.o object file (IBM compliler) or the datetime_g95.o file (g95 compiler). E-mail me both the source code and the location of the executable version of your program. ----------------------------- ARRAYS ------------------------------- Arrays are defined in F90 as in this example: real, dimension(100) :: a, b integer, dimension(50,2) :: index This defines a and b as 100 element arrays (a(1) to a(100)) and index as a 50x2 element array (index(1,1) to index(50,2)). Note that the default starting array index is 1 (not 0 as in C). Also note that array elements are written usng parentheses, not brackets as in C. Older Fortran programs would define the same arrays as follows: real a(100), b(100) integer index(50,2) This syntax will still work and is easier to read for short programs. The new standard has the advantage, however, of being able to easily define many arrays with identical dimensions. A nice feature of Fortran is that we can specify the lower and upper array boundaries explicitly in the declaration, e.g., real, dimension(-100:100) :: a, b integer, dimension(0:50, 2) :: index In this case there are 201 values for a and b (from a(-100) to a(100)) and 102 values for index (from index(0,1) to index(50,2)). Here is an example program that uses an array to compute prime numbers less than 100: program prime implicit none integer, parameter :: maxnum=100 integer :: i, j, prod(maxnum), max_i, max_j, nprime=0 do i = 1, maxnum prod(i) = 0 enddo max_i = floor(sqrt(real(maxnum))) do i = 2, max_i if (prod(i) == 0) then max_j = maxnum/i do j = 2, max_j prod(i*j) = 1 enddo end if end do do i = 2, maxnum if (prod(i) == 0) then nprime = nprime + 1 print "(i4)", i end if enddo print *, 'Number of primes found = ',nprime end program prime The method is sometimes called the sieve of Eratosthenes, named after a Greek mathematician from the 3rd century BC. We start with a list of numbers from 2 to 100 and consider them all possible primes. In the program, this list is the array prod. We initialize the array by setting its values to zero. The strategy is to then flag numbers that are not prime by setting the corresponding array values to one. We then start with the number 2 and eliminate all multiples of 2 up to the maximum value of 1000. We then move to 3 and eliminate (i.e., set prod to 1) all multiples of 3. We can skip 4 because 4 and all its multiples were already eliminated. We need check numbers, i, only up to sqrt(100) because larger factors would already have been eliminated. When we are finished, we simply print out all numbers that were not eliminated. In our case, this is all i such that prod(i) = 0. We count the number of primes using the counter variable nprime. A new aspect of this program is the defining of maxnum as a parameter: integer, parameter :: maxnum=100 This tells the program to set maxnum to 100 and that this value will never be changed during the program. It also allows the array dimension for prod (in the next line) to be set by maxnum. This makes it easy for us to change the size of our prime search by changing the value of maxnum without having to change anything else in the program. Note that in the statement max_i = floor(sqrt(real(maxnum))) it is necessary to convert maxnum to a real number before taking the square root. This is because the argument to the Fortran sqrt function must be real; if we had written sqrt(maxnum) we would have gotten an error during compilation. This program lists the prime numbers with one number per output line and thus will become cumbersome for larger values of maxnum. Let's modify the program to list 10 primes per line. To to this, we save the prime numbers in a separate array called pnum. Here is the code: program prime2 implicit none integer, parameter :: maxnum=1000 integer, dimension(maxnum) :: prod, pnum integer :: i, j, max_i, max_j, nprime=0 do i = 1, maxnum prod(i) = 0 enddo max_i = floor(sqrt(real(maxnum))) do i = 2, max_i if (prod(i) == 0) then max_j = maxnum/i do j = 2, max_j prod(i*j) = 1 enddo end if end do do i = 2, maxnum if (prod(i) == 0) then nprime = nprime + 1 pnum(nprime) = i end if enddo print *, 'Number of primes found = ',nprime print "(10i5)", (pnum(i), i=1,nprime) end program prime2 Note that we use the new F90 convention for setting up the arrays. The program is pretty self-explanatory, but we do introduce a new concept in the output line: print "(10i5)", (pnum(i), i=1,nprime) This is termed an implicit do loop and is very useful in I/O. Note that the format specifier 10i5 applies to the first 10 numbers and then is reused for the next 10, etc. The parentheses around (pnum, i=1,nprime) are required. More complicated expressions are also possible, e.g., print "(5(i4,i4,2x))", (i,pnum(i), i=1,nprime) will print the primes to the right of their count. Array values can be assigned one element at a time, or can be assigned in single statements as in this F90 example: program testarray implicit none integer, dimension(3) :: x = (/ 1, 2, 3 /), y = 1 print *, x x = (/ 15, 30, 40 /) print *, x print *, y y = 2 print *, y end program testarray and its output: 1 2 3 15 30 40 1 1 1 2 2 2 Note that when an array is set to a scalar, every array element is set to the value of the scalar, i.e., y =2 sets all y values to 2. The individual elements can be specified if they are listed between "(/" and "/)" (Warning: Don't put a space between the / and the parenthesis!!). The number of values must match the dimension of the vector. We will discuss this more later when we consider two-dimensional arrays in detail. So we see that we could have written: integer, dimension(maxnum) :: prod = 0, pnum in program prime rather than wasting the three lines we used later to set all elements in prod to zero. There are some cases when we would prefer not to advance to the next line following a print statement. Instead, we would like to have output continue on the same line. Here is an example of how this can be used to print 10 prime numbers per line without having to store the prime numbers in a separate array: program prime3 implicit none integer, parameter :: maxnum=1000 integer :: i, j, prod(maxnum)=0, max_i, max_j, nprime=0 max_i = floor(sqrt(real(maxnum))) do i = 2, max_i if (prod(i) == 0) then max_j = maxnum/i do j = 2, max_j prod(i*j) = 1 enddo end if end do do i = 2, maxnum if (prod(i) == 0) then nprime = nprime + 1 if (mod(nprime,10) /= 0) then write (*,"(i4)",advance='no') i else write (*,"(i4)") i end if end if enddo print * print *, 'Number of primes found = ',nprime end program prime3 The statement write (*,"(i4)") i works exactly the same as the print statement we used previously. ("print" only goes to the screen, "write" can go to the screen or to a file. The * in the above write statement means standard output, not free format). The statement write (*,"(i4)",advance='no') i does the same thing, except that it does not advance to the next line. The program checks to see if the prime count (nprime) is divisible by 10; if not then the output does not advance. Following the enddo, a print * statement is needed to be sure that the final output (Number of primes found) is on a separate line. NOTE: With the Suns, under both F77 and F90, the carriage return (advancement to next line) can also be suppressed simply by adding a $ to the format statement. Thus the above line could be replaced with print "(i4,$)", i and the result would be the same. However, I don't think this is standard Fortran so you may get into trouble at some point with this convention. ----- USING COMPILER OPTIONS TO CHECK FOR OVERFLOWING ARRAYS ------ Suppose when writing the prime number program, we made the mistake of writing: max_j = maxnum rather than max_j = maxnum/i When the program is run, at some point the product i*j in the line prod(i*j) will exceed the array dimensions of prod (maxnum). The program will then start overwriting adjacent memory locations and disaster will result. Most commonly the program will crash with the following message: rock% prime4 Segmentation Fault rock% Segmentation Faults or Bus Errors result from exceeding array boundaries or mismatched subroutine arguments. This is not a very helpful error message because it does not tell you where in the program the problem occurred. Many Fortran compilers provide a very useful option that provides better diagnostics. On the Suns, it is the -C option and can be invoked either by compiling prime2 with: f90 -C prime2.f90 -o prime2 or changing the Makefile to: %: %.f90 f90 -C $< -o $* (remember the tab before f90!) If you have compiled your program with the -C option, then the computer will check to see if your array indices exceed their boundaries and tell you where the problem occurs: ****** FORTRAN RUN-TIME SYSTEM ****** Subscript out of range. Location: line 15 column 20 of 'prime4.f90' Subscript number 1 has value 1002 in array 'PROD' Abort rock% For the G95 compiler on the Macs, use the -fbounds-check option, i.e., g95 -fbounds-check prime4.f90 -o prime2 or change the Makefile to: %: %.f90 g95 -fbounds-check $< -o $* (remember the tab before g95!) In this case, if you compile prime2.f90 with this option, the program will crash when executed and give the following message: At line 16 of file prime4.f90 Traceback: not available, compile with -ftrace=frame or -ftrace=full Fortran runtime error: Array element out of bounds: 1002 in (1:1000), dim=1 These messages are so helpful, and will save you so much time in debugging your programs, that I recommend that you ALWAYS use these options when you are first getting your programs to work. It does slow the code somewhat, so for programs where run time is a factor, you should then recompile without the -C (or -fbounds-check) option for the final working version. (in these cases, you may also want to experiment with the -O compiler to make your code run faster, see below). IBM compiler: Under IBM F90, the prime4 program above will give "Bus error" during run time. The -C option can also be used with the IBM F90 compiler so that the program will terminate during run time if an array goes out of bounds. In this case the prime4 program above will terminate with "Trace/BPT trap" during run time. Unfortunately, it does not give you the details (line number, etc.) the way the Sun F90 compiler does. However, it appears from the documentation (http://www.absoft.com/Products/Compilers/Fortran/Macintosh/XLF/xlf/pgs/ug31.htm) that it is possible to get more information if one understand how SIGTRAP works. This looks complicated to me, but it would be great if someone figured this out. Info seems to be at: http://www.absoft.com/Products/Compilers/Fortran/Macintosh/XLF/xlf/pgs/ug32.htm#HDRSIGTRAP for the -qsigtrap option, and http://www.absoft.com/Products/Compilers/Fortran/Macintosh/XLF/xlf/pgs/ug39.htm#HDRHANDLER for the "Installling an Exception Hander" section. It appears that this may involving putting some additional lines into the Fortran source code. Even without more details, it is still helpful to get the "Trace/BPT" error so I recommend that you use the -C option until you have the final working version of your code. ASSIGN F11 Fortran 90 includes a built in random number subroutine (called "random_number"). Here is an example program that prints out 20 random values from between 0 and 1. program listrand integer :: i real :: xran do i = 1, 20 call random_number(xran) print *, xran enddo end program listrand Write a simple F90 program to generate 10000 random numbers between 0 and 1. Test the randomness by counting the number of these that fall in each of 10 evenly spaced bins between 0 and 1 (i.e., 0 to 0.1, 0.1 to 0.2, etc.). Print these counts to the screen. Here is an example of what your output might look like: 0 1009 1 1048 2 1001 3 1038 4 1008 5 959 6 993 7 925 8 1017 9 1002 HINT: Set up an integer array with 10 elements and initialize all the values to zero. For each random number generated, then add one to the appropriate array element to keep track of how many of the numbers are in that bin. NOTE: You will get the same random numbers each time with random_number. In principle, you get change this by using random_seed in F90 but this is complicated to code (see http://www.atmos.umd.edu/~miyoshi/memorandum/memo_random.html). A better choice for serious work are the random number generators in Numerical Recipes. ---------------- ARRAYS AS SUBROUTINE ARGUMENTS ------------- Suppose we have defined an array x as: integer :: x(100) If we use x in the argument list for a subroutine, e.g., call SUMTOT(x,n,y) we include x without any parentheses. The subroutine argument list must also include a matching array of the same size. All array values are passed to and from the subroutine. We will discuss later how to write subroutines that will work for arrays of differing size in the calling program. --------------------- CHARACTER STRINGS ------------------------ A character string may be declared in F90 as follows: character (len = 20) :: name declares the variable "name" to be a character string of length 20. One can also define an array of character strings: character (len = 20), dimension(100) :: name_array The old F77 ways to define these strings will still work in F90: character name*20, name_array(100)*20 OR character*20 name, name_array(100) A string variable can be assigned as, e.g., name = "Bill Clinton" The quotes (apostrophes will also work) are required. If name in this case was declared as 20 characters long, then trailing blanks are added so name actually is "Bill Clinton " If, on the other hand, name was declared as 10 characters long, then name would contain "Bill Clint" as the extra letters would be cut off (no error results). The length of a character string may be obtained with the len function: len("Bill") will give 4 len(" ") will give 2 len(name) will give 20 if name was declared as 20 characters long, regardless of any trailing blanks Substrings may be obtained as follow for name = "Bill Clinton": name(1:4) is "Bill" name(6:12) is "Clinton" name(6:6) is "C" Thus, we could change the last name as follows: name = "Bill Clinton" name(6:12) = "Jones " Note that the trailing blanks are necessary to overwrite the "on" The following moves the string one character to the right: name(2:13) = name(1:12) This was not allowed in F77 (because of the overlap) but is OK in F90. Often one wants to know the length of a string without any trailing blanks. This can be done with the built in function len_trim. Another built in function is called "trim" which gives the string without any trailing blanks. To find the position of a substring within a string, the built in funtion index can be used: index("Clinton", "in") returns 3 index("Clinton", "on") returns 6 index("Clinton", "no") returns 0 (indicates string not found) To append one string onto the end of another (concatenation), the // operator can be used: text1 = "Bill Clinton" // " was elected in 1992." sets text1 to the complete sentence (text1 must have been already declared of sufficient length). Note the blank before "was" is necessary for there to be a space between the word. When trailing blanks are present, the "trim" function is very useful. For example, if name is declared as 20 characters long, then name = "Bill Clinton" text1 = name // " was elected in 1992" sets text1 to "Bill Clinton was elected in 1992" while text1 = trim(name) // " was elected in 1992" will produce the correct result. Here is an example program that illustrates how to input a string: program vote implicit none character (len = 80) :: name print *, 'Who did you vote for in 1992?' read "(a)", name if (index(name,'ill') /= 0 .or.index(name,'lin') /= 0) then print *, "Then you are likely a Democrat." else if (index(name,'eor') /= 0 .or.index(name,'ush') /= 0) then print *,"Then you are likely a Republican." else print *,"Then you are likely an independent voter." end if end program vote Note the format for the read statement: read "(a)", name does not need to use a80 although that would also work. The free format read statement: read *, name is not recommended in this case because it will only read until the first blank. Thus Bill Clinton would be read in as "Bill" (unless the full name were enclosed in quotes, but who wants to make the user do that?). ASSIGN F12 Write a simple version of the famous doctor program that simulates a psychologist talking to a patient. A humorous fictional account of such a program is contained in "Small World," the David Lodge satire of academia. A current example of this type of program can be found at alicebot.org (click on "Chat with Alice"). Begin by asking the patient: "How are you?" Then use the index function to input a line from the user/patient. Search this line for various key words that you can then use to guide the computer's response. Design the program so that it will continue this "conversation" between doctor and patient. Some suggestions: If the input string contains a "?", then respond: "Let me ask the questions. Why is it important to you?" If the string contains "!", then respond: "You seem pretty upset. What is really bothering you?" If the string contains "mother" then respond: "Tell me more about your mother." If the string contains any swear words, respond: "There is no need for that kind of language." If you identify no key words on your list, then respond with something generic like: "Go on." or "Tell me more about it." Be creative but don't spend an infinite amount of time on this assignment. Your program will be more realistic if you use random numbers (see above) to randomly select from a series of possible responses so that the computer does not keep saying the same thing. ------------------------- I/O WITH FILES -------------------------- So far, all off our examples have involved input from the keyboard and output to the screen. Now let us see how to open files, read from them, and write data to them. Here is a program that reads pairs of numbers from an input file, computes their product, and outputs the original numbers and their product to an output file. program fileinout implicit none character (len=100) :: infile, outfile real :: x, y, z integer :: ios print *, 'Enter input file name' read "(a)", infile open (11, file=infile, status='old') print *, 'Enter output file name' read "(a)", outfile open (12, file=outfile) do read (11,*, iostat = ios) x, y if (ios < 0) exit z = x * y write (12,*) x, y, z enddo close (11) close (12) end program fileinout We open the input file with the statement: open (11, file=infile, status='old') Files must be assigned "unit numbers" for I/O. This statement opens the file with unit 11. On most systems, unit 5 is defined as the default standard (keyboard) input and unit 6 is defined as the default standard (monitor) output, so these numbers should always be avoided as unit numbers. I have gotten into the habit of using units 11, 12 and 13 for most of my I/O. file=infile is what provides the filename. Note that we could "hardwire" a file name here by writing: open (11, file="file1", status='old') The status='old' is optional, but I recommend always using it for opening files which should already exist. When you say status='old' and the file does not exist, then the program will terminate and print an error condition. If you have not set status='old' then the program will create a new file and when you try to read from it, you will reach the end of the file immediately and the program will terminate without any obvious errors (except that you will now have zero length files in your directory!). Another value for status= is: status='new' file must not already exist (useful to avoid overwriting existing files) However, in this example we choose to simply write: open (12,file=outfile) so that we can overwrite existing file names. This assigns unit 12 to the output file (the input and output file unit numbers must be different!). We read data from the input file with the statement: read (11, *, iostat = ios) x, y This does a free-format read from unit 11 (for a fixed format, the * is replaced with a format specifier such as "(i2,2f6.2)" etc.). The only complication is that we need some way to recognize when we have reached the end of the file (EOF). The "iostat = ios" specifier does this by assigning the local variable ios (which we must declare as an integer) that is zero for a normal read, positive if there is an error during the read, and negative if the end of the file is reached. Note that we are free to choose any name for ios that we want as long as it is declared as an integer. Next, we use if (ios < 0) exit to exit the otherwise infinite do loop when the EOF is reached. We should note at this point that older Fortran programs do not use this method of testing for the end of the file. Rather, they will say something like: read (11, *, end = 123) x, y where 123 is a statement label that the program will branch to when the end of file is reached. This method still works in F90 but you lose style points for having to use a statement label. We write the original numbers and their product using a free format write to unit 12: write (12,*) x, y, z Note that x and y will probably not be in the same format as in the input file. Finally, after the EOF but before we end the program it is good practice to close the files: close (11) close (12) although the files will be closed automatically anyway when the program ends. ASSIGNMENT F13 Write a program that reads from a text file containing 5 numbers per line. Compute the mean of the 5 numbers and then subtract the mean from each of the 5 original numbers. Output the 5 "demeaned" numbers to an output file. Allow the user to specify the input and output file names. For example, if the input file contains: 10 20 30 40 50 1 2 3 4 5 2 4 6 8 10 5 5 5 5 5 your program should write to the output file: -20.000 -10.000 0.000 10.000 20.000 -2.000 -1.000 0.000 1.000 2.000 -4.000 -2.000 0.000 2.000 4.000 0.000 0.000 0.000 0.000 0.000 Hint: Use the "fileinout.f90" program listed above as a template for the file handling part of your code. --------------------- TWO-DIMENSIONAL ARRAYS -------------------- Two-dimensional arrays in Fortran may be defined as follows: integer a(2,3) !F77 convention integer, dimension(2,3) :: a !F90 convention This defines a matrix in which the first indices will vary from 1 to 2, and the second from 1 to 3 (recall that the default starting array index is 1 in Fortran). Here is a test program that illustrates how the numbers in a can be specified and how they are stored in memory: program testmatrix implicit none integer, dimension(2,3) :: a a(1, 1:3) = (/ 1, 2, 3 /) a(2, 1:3) = (/ 4, 5, 6 /) print *, a end program testmatrix The output of this program is: 1 4 2 5 3 6 First, note how the array values are specified. (/ 1, 2, 3 /) is an example of what is called an "array constructor" and is a sequence of scalar values along one array dimension only. We equate this to the desired part of the a matrix by writing, for example: a(1, 1:3) = (/ 1, 2, 3 /) thus setting a(1,1)=1, a(1,2)=2, and a(1,3)=3. Finally note how the array is printed. The print *, a statement will dump the values of a in the same order that they are stored in memory. In doing this, Fortran varies the first array element first, then the second, etc. In this case, the result is: a(1,1), a(2,1), a(1,2), a(2,2), a(1,3), a(2,3) Similarly a three-order array b(2,2,2) would be stored as: b(1,1,1) b(2,1,1) b(1,2,1) b(2,2,1) b(1,1,2) b(2,1,2) b(1,2,2) b(2,2,2) This is important to remember when passing arrays to and from subroutines for cases when the array in the main program has more dimensions than the array in the subroutine. For example, suppose we wished to store 100 seismograms of 1000 points each in the main program. We have a filtering subroutine that will act on a single time series vector. If we dimension the array in the main program as b(1000,100) then we can call the subroutine using a statement like: call FILTER(b(1,i),npts) where i is the seismogram number. The filter can then act on b(1:npts, i) and would be of the form: subroutine FILTER(c, n) real, dimension(n) :: c . . This is a very common type of construction. Note that this would not work if the array was dimensioned as b(100,1000) because then each 1000-point time series would not be continuous in memory. Here is an example program to multiply two matrices: program matmult implicit none real, dimension(3,3) :: a, b, c integer :: i, j, k a(1,1:3) = (/ -5.1, 3.8, 4.2 /) a(2,1:3) = (/ 9.7, 1.3, -1.3 /) a(3,1:3) = (/ -8.0, -7.3, 2.2 /) b(1,1:3) = (/ 9.4, -6.2, 0.5 /) b(2,1:3) = (/ -5.1, 3.3, -2.2 /) b(3,1:3) = (/ -1.1, -1.8, 3.0 /) do i = 1, 3 do j = 1, 3 c(i,j) = 0.0 do k = 1, 3 c(i,j) = c(i,j) + a(i,k) * b(k,j) enddo enddo enddo print *, "Matrix a follows" call PRINTMAT(a) print *, "Matrix b follows" call PRINTMAT(b) print *, "Matrix c = a*b follows" call PRINTMAT(c) contains subroutine PRINTMAT(x) real, dimension(3,3) :: x do i = 1, 3 print "(3f8.3)", (x(i,j), j=1,3) enddo end subroutine PRINTMAT end program matmult Note the use of the internal subroutine PRINTMAT to display the results. This example performs the matrix multiplication exactly as one would in F77; it does not take advantage of some of the new array capabilities of F90. In fact, there is a built in function, MATMUL, in F90 that performs matrix multiplication. My plan is to have more to say about these special F90 capabilities later in the class. ASSIGNMENT F14 Modify the matmult program so that the matrix multiplication is done in a subroutine. Add another subroutine that computes the transpose of a square matrix. Compute the transpose of c and add this to the output listing. ---------------------- ARRAYS OF STRINGS ----------------------- As we have discussed, a string is just an array of characters. Thus, a 1-D array of strings will be a 2-D character array. Here is an example of this: program stringarray implicit none character (len = 80), dimension(10) :: name name(1) = "Bill Clinton" name(2) = "George Bush" print *, "Our last president was ", trim(name(1)) print *, "Our current president is ", trim(name(2)) end program stringarray The program can store up to 10 name of 80 characters each. Note the use of the "trim" function when printing the names. This is very handle to cut off the trailing blanks and limit unnecessary extra spaces (and sometimes lines) in the output. --------- A MORE COMPLICATED EXAMPLE OF PROCESSING A DATA FILE ---------- Suppose we have a data file that contains a series of measurements that has the following form: (datafile1) 980204 22:03:34.42 10.8 10.2 9.8 11.6 980205 08:45:22.20 5.5 7.2 6.6 980205 20:13:42.88 15.0 13.8 12.9 14.2 etc. The structure is that there are date/time identifiers and then a series of measurements made at that time. The problem in reading this file is that we don't know how many measurement lines will follow each date/time line so we can't simply perform a formatted read on each line from the file. Let us assume that we want to compute the mean of these measurements and write a new file which stores this information on the same line as the date.time. In this case, our desired output file will look like: (datafile2) 980204 22:03:34.42 10.60 4 980205 08:45:22.20 6.43 3 980205 20:13:42.88 13.97 4 where the final column gives the number of measurements. Here is how such a program could be written: program procfile implicit none character (len=100) :: infile, outfile character (len=20) :: linebuf, event integer :: ios, nevent=0, ndata=0 real :: x, sum=0.0, xavg print *, 'Enter input file name' read "(a)", infile open (11, file=infile, status='old') print *, 'Enter output file name' read "(a)", outfile open (12, file=outfile) do read (11,'(a)', iostat = ios) linebuf if (ios < 0) exit if (linebuf(1:1) /= ' ') then ! no init blank, must be event nevent = nevent + 1 if (nevent /= 1) then !need to output previous event info call WRITE_EVENT ndata = 0 sum = 0 end if event = linebuf else ! has starting blank, must be measurement read (linebuf,*) x sum = sum + x ndata = ndata + 1 end if enddo call WRITE_EVENT close (11) close (12) contains subroutine WRITE_EVENT if (ndata /= 0) then xavg = sum/real(ndata) else xavg = 0. end if write (12,"(a20, f7.2, i5)") event, xavg, ndata end subroutine WRITE_EVENT end program procfile The beginning of the program where the files are opened is the same as in fileinout.f90 (our earlier example of file I/O). Using an infinite do loop, we read lines from the input file using: read (11,'(a)', iostat = ios) linebuf Because we don't know the format of each line before we read it, we read into a character string (linebuf) that will serve to temporary store the contents of the line (a line "buffer" as they say in the computer world). We use iostat=ios to set ios to a flag that can be used to identify when the end of the file occurs so that we can exit the do loop. Next, we check the contents of the first character in linebuf. If is not a blank, then we know we have an event line. In this case we increment the nevent counter variable by one. If we have not just read the first event, then we have read a new event, thus terminating the data input from the previous event. In this case we need to write the processed event info to the output file (we can't do this until this point because we don't know how many data points there will be) and reset sum and ndata to zero. Finally we set event = linebuf to store the event line until we need to output it. If, on the other hand, the first character in linebuf is blank, then we have a data line. We then read the value of x out of linebuf: read (linebuf,*) x This is called an "internal read" and is the first time that we have shown this. It allows the program to do a formatted read from a string, rather than from an external file. We do a free format read (*), but one could also use a format specifier. After reading the measurement value into x, we suitably increment a running sum of x and the number of measurements. When we reach the end of the input file, we have to output the results of the final set of measurements. To avoid repeating a block of code for this operation, we use an F90 internal subroutine (WRITE_EVENT). This is uses twice, inside the loop as we output information from the last event before going on to the next one, and outside the loop where we output the information from the last event. Note that we must check to see if there are no measurements for an event line, in which case we do not divide by zero but simply set xavg to zero. There should be no possibility of mistaking this for a real measurement in the output file because the number of measurements in this case will also be zero. ASSIGNMENT F15 Obtain the data file, scsn.phase.dat, from mahi.ucsd.edu/pub/239. This is a list of earthquakes and arrival time data from the Southern California Seismic Network (SCSN) for the first two days of August 2000. (Files like this are readily obtainable from the SCEC Data Center at http://www.scecdc.scec.org if you ever want to obtain more of these data. The format is called the SCEC_DC phase format.) Here is what part of the file looks like: 2000/08/01 14:44:12.6 L 1.1 c 35.962 -117.693 5.2 A 9159024 8 51 0 0 WRC VHZ 911 P IU1 4.41 1.26 WCS VHZ 885 P E 3 9.98 1.85 WVP VHZ 925 P IU1 11.54 2.25 TOW VHZ 813 P E 2 18.28 3.78 CLC VHZ 202 P E 2 18.49 3.37 WMF VHZ 906 P E 3 22.86 3.94 WMF VHZ 906 S E 2 22.86 7.08 WNM VHZ 908 P E 3 23.74 4.23 2000/08/01 14:45:36.8 L 1.3 h 34.806 -116.271 7.3 A 9159025 7 79 0 0 CDY VHZ 184 P IU0 6.66 1.73 RAG VHZ 1034 P IU1 17.58 3.19 RAG VHZ 1034 S E 2 17.58 5.77 RMM VHZ 647 P ID1 37.23 6.33 TPM VHZ 1033 P IU1 52.29 8.88 GRP VHZ 355 P IU1 61.21 10.14 0299 C 65535 P IU1 11395.79 7.39 2000/08/01 15:04:17.9 L 1.4 h 34.810 -116.268 7.1 A 9159030 19 83 0 0 CDY VHZ 184 P IU0 6.76 1.73 RAG VHZ 1034 P IU1 17.76 3.16 etc. There is a line for each earthquake with the date and time (UT not local time). For the first line in the above example, the additional bits include: L = local event 1.1 = magnitude c = how magnitude was computed (coda in this case) 35.962 = lat -117.693 = lon 5.2 = depth (km) A = SCSN assigned quality (A=best) 9159024 = cusp id number 8 = number of lines of phase arrival time info The station info lines include: WRC = station name VHZ = station channel 911 = ? P = phase name (normally P or S) IU1 = phase pick info 4.41 = station-event distance (km) (=X for assignment) 1.26 = travel time (s from event origin time) (=T for assignment) Your task is to write a F90 program to read this file and write to an output file a series of (X,T) points with (one X-T pair per line) for all P arrivals for all events between 3 and 8 km depth. Then make an X-Y plot of the T-X points using PLOTXY. Finally, write a UNIX script that will run the program, run PLOTXY and display the plot on the screen, all by simply typing the name of the script. E-mail me a copy of the F90 program source code and the UNIX script. The good news about this data format is that it does include a number in the event line that lists the number of phase lines that will follow. The bad news is that sometimes the last 1 or 2 phase lines are "garbage" in that they have numbers instead of station names and are some kind of calibration info. You will have to figure out how to recognize and discard these lines in your program. Also the phase lines begin with a tab character rather than a long series of blanks; this may complicate things somewhat. Hint: Don't try to write your entire program all at once. First, get a version working that simply reads the input file and prints out the parts of each line that you will need (i.e., the depth and number of phase lines for the event line, the phase name and the X and P values for the phase lines). Once you have this working correctly, then go on to add the part to output the (X,T) points. -------------- SORTING EXAMPLE USING NR UTILITIES --------------- The book Numerical Recipes contains a large number of useful subroutines (there are both Fortran and C versions) that are fully explained in the text. You can find F77 source code (1st edition of NR) for these routines in: /net/rock/shearer/SUBS/NUMRECIP To illustrate the use of the Numerical Recipes routines, here is an example program that generates a list of random numbers, sorts them using the NR function piksrt, and then prints out the sorted list. program sortrand implicit none integer, parameter :: NPTS=100 integer :: i real, dimension(npts) :: xran do i = 1, NPTS call random_number(xran(i)) enddo call PIKSRT(NPTS,xran) print "(10f7.4)", (xran(i), i=1,npts) end program sortrand The program sets the xran array to random numbers. It then calls the NR subroutine PIKSRT to sort the numbers before printing them out. We must either include the PIKSRT source code or link with a PIKSRT object code when we compile this program. Here is the PIKSRT source code: SUBROUTINE PIKSRT(N,ARR) DIMENSION ARR(N) DO 12 J=2,N A=ARR(J) DO 11 I=J-1,1,-1 IF(ARR(I).LE.A)GO TO 10 ARR(I+1)=ARR(I) 11 CONTINUE I=0 10 ARR(I+1)=A 12 CONTINUE RETURN END This ugly code is typical of old Fortran programs. It has all capital letters and has line numbers, "go to" statements, and uses "continue" statements in the do loops. Actually, it's better than many examples, because the do loops are indented. The important thing for us is that the code works--we don't want to bother rewriting it and possibly breaking it. In general, F77 and older code will not compile under F90 because of old-style comment lines ("C" in column 1) and old-style line continuation flags (non-& character in column 6). This example, however, does not have these problems and is fully F90 compatible. Thus, we could just append the source code onto the end of our program. Alternatively, we could compile the subroutine under F90 as follows: f90 -c piksrt.f -o piksrt_f90.o which will create a piksrt_f90.o file. The _f90 is my own convention to label f77 source code that has been compiled under f90 so that I don't try to link piksrt_f90.o with a f77 program. Unfortunately, f77 and f90 object files are not always compatible (they are in this case, but it is dangerous to assume that they always will be, so it is safest to separately compile the source code in each case). f90 sortrand.f90 piksrt_f90.o -o sortrand or by using a Makefile as discussed earlier. This would be our best choice if the subroutine source code is lengthy and/or is non-F90 compatible. ASSIGNMENT F16 Modify the sortrand program to return the median of a user-specified number of random numbers. Do the median computation in the form of a subroutine that returns the median value, but does NOT resort the input array in the process. The random number generation should be in your main program, not within the MEDIAN subroutine which should be a general purpose routine that simply inputs an array of numbers and returns the median value. Then whenever you need to compute a median of some numbers within a program, you can simply call the same subroutine, e.g., call MEDIAN(x, n, xmed) !x=array, n=# of points, xmed=median Hint: Within your median-computing subroutine, copy the input array to another array before calling PIKSRT. Make sure that your program gives the correct result for both even and odd numbers of points. -------------- EXAMPLE OF SAVING VALUES IN A SUBROUTINE -------------- A common convention in geophysics is to name an instrument site with a short ascii string. In seismology, for example, station names are usually designated with 3 to 4 character names. Southern California GPS sites are usually identified with 4 character names. Often data products are labeled by the station name, without including information about the station, such as its location. Data processing programs will often require this information. Thus, there is frequently a need for a computer routine that will retrieve information about a station, given the station name. One could hardwire this information into the program with a series of if statements, but this would become very cumbersome for large numbers of stations and would limit the flexibility and portability of the code. A better approach is to save the station information in a file and have the computer read from the file to access the information when necessary. However, file I/O is relatively slow so we won't want to open, read, and close the file every time we need to determine a station location. It would be much faster to read the station information file once and then save the information during the time that the program is running. For portability, ideally all of this overhead would be performed in a function that we could use in many programs without having to worry about the format of the station file or the size of the array that are required to store it. In order to do this, we need to retain the values of certain variables within the function for subsequent calls. Here is an example that will return station coordinates for some of the GSN (Global Seismic Network) stations. program testgetstat implicit none character (len = 4) :: stname real :: slat, slon, selev do print *, 'Enter station name (stop to stop)' read (*,'(a)') stname if (stname == 'stop') exit call GET_STAT(stname,slat,slon,selev) print *, 'slat,slon,selev = ',slat,slon,selev enddo end program testgetstat ! subroutine GET_STAT gets the lat/lon/elev of a station ! with a given name from the GSN station list ! ! Inputs: snam = station name (a4) ! Returns: flat = station latitude ! flon = station longitude ! felev = station elevation (m) ! NOTES: Station list must be alphabetized ! If station name is not found, then returned ! variables are set to -999. ! Program reads and saves station info on first call ! subroutine GET_STAT(snam,flat,flon,felev) implicit none integer, parameter :: NMAX=5000 character (len = 4) :: snam character (len = 4), dimension(NMAX) :: stname real, dimension(NMAX) :: slat, slon, selev real :: flat, flon, felev integer :: i, i1, i2, nsta, it logical :: firstcall = .true. save firstcall, stname, slat, slon, selev, nsta if (firstcall) then firstcall = .false. open (11, file='stlist.gsn', status='old') print *, 'Reading station file: stlist.gsn' do i=1,NMAX read (11,7,end=12) stname(i), slat(i), slon(i), selev(i) 7 format (a4, f10.5, f11.5, f5.0) enddo print *,'***Warning: number of stations in file may exceed ',NMAX i = NMAX + 1 12 nsta = i-1 close (11) print *,'Number of stations read = ',nsta end if i1 = 1 i2 = nsta do it = 1, 15 i = (i1+i2)/2 if (snam == stname(i)) then flat = slat(i) flon = slon(i) felev = selev(i) return else if (snam < stname(i)) then i2 = i-1 else i1 = i+1 end if enddo print *, '***station not found ', snam flat = -999. flon = -999. felev = -999. end subroutine GET_STAT The main program is a short driver program that enables us to test that the GET_STAT subroutine is working properly. It is always a good idea to use little programs like this to test functions that are being developed, BEFORE using the functions in larger programs. For efficiency, the subroutine reads the station file only the first time the subroutine is called. It then saves the information in memory to use for subsequent calls. Normally, the values of variables in a subroutine are not retained between calls. The "save" command is used to specify those variables that are to be saved between calls. This is our first example of a logical variable: logical :: firstcall = .true. firstcall can have two possible values: .true. or .false. Here we set it to .true. the first time (and only the first time) that the subroutine is called. The save statement specifies which variable values are to be saved between subroutine calls: save firstcall, stname, slat, slon, selev, nsta To read the station file upon the first call to the subroutine, we write if (firstcall) then Notice that a logical variable can be the sole argument for an if test. We then set firstcall = .false. so that this block of code will not be excuting upon subsequent calls to the subroutine. The station file has the form: AAE 9.02920 38.76560 2442 AAK 42.63900 74.49400 1645 ABKT 37.93040 58.11890 678 ADK 51.88370 -176.68440 116 AFI -13.90930 -171.77730 706 ALE 82.50330 -62.35000 60 ALQ 34.94620 -106.45670 1840 ANMO 34.94620 -106.45670 1840 ANTO 39.86890 32.79359 883 AQU 42.35389 13.40500 720 . . In routines like this, I like to print out information about the file that is being read during the first call. This reminds the user of the program about what is being done and it helpful in case there is a problem or error later in the program. For example, perhaps we have a faulty station file that has only 1 station. Or perhaps the firstcall flag is not working properly, in which case these output lines will result each time we call the routine, not just the first time. For speed, the subroutine does not simply loop through the station list until it finds a match. Instead, it exploits that fact that the stations are in alphabetical order. i1 and i2 are indices that are designed to bracket the target station in the list. Initially, they are set to 1 and nsta, the first and last station indices. Next, i is set to a point halfway between i1 and i2. The station in the station list at i is compared to the target station. If they are the same, then we can set the return variables. If the stlist station is greater than the target station (later in the alphabet), then i2 is set to i-1. If the stlist station is less than the target station, then i1 is set to i+1. We iterate using this procedure 15 times to narrow in on the station name. (2^15 must be greater than NMAX to make sure we do enough iterations) Note that else if (snam < stname(i)) then is a valid check to see if string variables are in alphabetical order. If we don't find the target stname in the station list, then we print out a warning message. We also set the station coordinates to numbers that are unlikely to be confused with real station locations in case the user of the function does not notice the warning message (or chooses to ignore it!). ASSIGNMENT F17 Modify GET_STAT to find the nearest station to a given (lat,lon) point and return the station name and its coordinates. You can get the stlist.gsn file via anonymous ftp at mahi.ucsd.edu/pub/239. Name your new function NEAR_STAT and also include a testnearstat driver program to test the operation of the function. To find the nearest station, you will need to be able to compute the distance between two points on a sphere. Use the SPH_AZI from the notes or the SPH_DIST subroutine from ~shearer/SUBS/sphere_subs.o Don't bother to do the separate calculation for small values of del. The NEAR_STAT suboutine should be of the form: subroutine NEAR_STAT(plat, plon, stname, slat, slon, sdep) where plat and plon are the coordinates going into the subroutine and stname, slat, slon and sdep are passed back from the function to the main program. ------------------------- COMPLEX NUMBERS -------------------------- A nice feature of Fortran is that complex numbers are a standard feature. Here is a program to show how they can be declared and used: program testcomplex implicit none complex :: a, b, c print *, 'Enter first complex number' read *, a print *, 'Enter second complex number' read *, b c = a*b print *, 'Product = ', c print *, 'abs of product = ', abs(c) print *, 'sqrt of product = ', sqrt(c) end program testcomplex Here is an example of running the program: rock% testcomplex Enter first complex number (1, 1) Enter second complex number (-.5, .1) Product = (-0.6,-0.4) abs of product = 0.7211103 sqrt of product = (0.2460795,-0.81274545) Complex numbers are represented as a pair of numbers in parenthesis separated by a comma (first number is the real part, second number is the imaginary part). This is the only format that will work for the free format read. If you try to enter 1, 1 or 1 1 or even (1 1) you will get an error. Notice that in F90 you don't need to use special forms for the functions abs and sqrt (e.g., cabs and csqrt as you will see sometimes in older code). There are conversion functions to build complex numbers from two real numbers and vice versa, as shown in this example code: program testcomplex2 implicit none complex :: a real :: ar, ai print *, 'Enter real part' read *, ar print *, 'Enter imaginary part' read *, ai a = cmplx(ar, ai) print *, 'a = ',a a = exp(a) print *, 'Exp(a) = ', a ar = real(a) ai = aimag(a) print *, 'Real part = ', ar print *, 'Imag part = ', ai end program testcomplex2 and its output rock% testcomplex2 Enter real part 1 Enter imaginary part 2 a = (1.0,2.0) Exp(a) = (-1.1312044,2.4717266) Real part = -1.1312044 Imag part = 2.4717266 ASSIGNMENT F18 Write a program to test the built in complex exponential function against the formula e^z = e^(x+iy) = e^x e^(iy) = (e^x cos y) + i (e^x sin y) which you will have to adapt for your program. Include a driver program to test your function for selected input values. Make sure that your code gets the correct answer for these examples: exp( 1.1 + 2.3i) = -2.002 + 2.240i exp( 0.0 + 1.2i) = 0.362 + 0.932i exp(-0.5 + 0.0i) = 0.607 + 0.000i Can you find any differences between the built in exp function are the results obtained with the formula? -------------------- F90 ARRAY ARITHMETIC ------------------- F90 allows many operations on vectors and matrices to be performed in single statements without the necessity of writing do loops over the array indices (as was necessary in F77). Here is an example program that demonstrates some of these operations on vectors. program vectormath implicit none real, dimension(5) :: a = (/ 1.0, 2.0, 3.0, 4.0, 5.0 /), b = 2.0, c real :: x print *, 'a = ', a c = a + 1 print *, 'a + 1 = ', c c = 2 * a print *, '2 * a = ', c c = a * a print *, 'a * a = ', c c = sqrt(a) print *, 'sqrt(a) = ', c c = sin(a) print *, 'sin(a) = ', c c = exp(a) print *, 'exp(a) = ', c print *, 'b = ', b c = a + b print *, 'a + b = ', c c = a * b print *, 'a * b = ', c x = sum(a) print *, 'sum(a) = ', x c = a c(4:5) = 0.0 print *, 'a with two zeros on end = ', c x = dot_product(a, b) print *, 'a dot b = ', x x = sum( (a - sum(a)/5. )**2) print *, 'sum of squares of difference from mean = ', x end program vectormath Note that even fairly complicated expressions are possible provided the programmer keeps track of what is a vector and what is a scalar. Operations are also possible on matrices. Here are some examples: program matrixmath implicit none real, dimension(2, 3) :: a23, b23, c23 real, dimension(3, 2) :: a32, b32, c32 real, dimension(2,2) :: a22, b22, c22 integer, dimension(2) :: loc2 integer :: i, j, k a23(1,1:3) = (/ -5.1, 3.8, 4.2 /) a23(2,1:3) = (/ 9.7, 1.3, -1.3 /) print *, 'Matrix a23 follows' call PRINTMAT(a23, 2, 3) b32(1:3, 1) = (/ 9.4, -6.2, 0.5 /) b32(1:3, 2) = (/ -5.1, 3.3, -2.2 /) print *, 'Matrix b32 follows' call PRINTMAT(b32, 3, 2) c22 = matmul(a23, b32) !this works but matmul(b32, a23) does not print *, "Matrix c22 = matmul(a,b) follows" call PRINTMAT(c22, 2, 2) print *, 'maxval(a23) = ', maxval(a23) print *, 'maxloc(a23) = ', maxloc(a23) loc2 = maxloc(a23) print *, 'loc2 = ', loc2 print *, 'a23(loc2(1), loc2(2)) = ', a23(loc2(1), loc2(2)) b23 = a23 + transpose(b32) print *, 'Matrix b23 = a32 + transpose(b32) follows' call PRINTMAT(b23, 2, 3) print *, 'size(a23) = ', size(a23) print *, 'size(a23(1,:)) = ', size(a23(1,:)) contains subroutine PRINTMAT(x, m, n) real, dimension(m, n) :: x integer :: m, n do i = 1, m print "(10f8.3)", (x(i,j), j=1,n) enddo end subroutine PRINTMAT end program matrixmath This demonstrates the intrinsic functions matmul, maxval, maxloc, transpose, and size. Note that matmul(a,b) is true matrix multiplication, not the multiplication of the individual array elements as occurs from writing simply a*b. The maxloc function returns the location of the maximum value in the array. Note that the statement loc2 = maxloc(a23) works only if loc2 is defined as an integer array with dimension(2). NOTE: If maxloc is applied to a one-dimensional array, for example to find the index of the maximum value of a vector bvec: loc = maxloc(a) the loc must be defined as an integer array with dimension(1). You will get an error message during compilation if loc is defined simply as an integer. There are corresponding functions, minval and minloc, that determine the minimum value and location within an array. It would be interesting to test whether programs using the intrinsic array functions in F90 run any faster than those written using do loops. ------------------------ ALLOCATABLE ARRAYS ---------------------- An awkward aspect of older Fortran programs is the need to declare arrays to the maximum possible size that could ever be needed when running the program. This memory must be set aside even when it is not needed--a wasteful practice. F90 avoids this by allowing memory to be allocated and deallocated on the fly. Here is an example program: program setarray implicit none real, dimension(:), allocatable :: x real, dimension(:, :), allocatable :: y, yy real :: xsum integer :: n, i, j, m print *, 'Enter number of points' read *, n allocate (x(n)) do i = 1, n call random_number(x(i)) enddo xsum = sum(x) print *, 'sum = ', xsum deallocate(x) !this frees up memory print *, 'Enter m, n (# rows, # columns)' read *, m, n allocate (y(m,n)) allocate (yy(n,n)) do i = 1, m print *, 'Enter row # ', i read *, (y(i, j), j = 1, n) enddo yy = matmul(y, transpose(y)) print *, 'y * yt = ' call PRINTMAT(yy, m, m) deallocate(y) deallocate(yy) contains subroutine PRINTMAT(x, m, n) real, dimension(m, n) :: x integer :: m, n do i = 1, m print "(10f8.3)", (x(i,j), j=1,n) enddo end subroutine PRINTMAT end program setarray -------------------- STRUCTURES -------------------- It is often convenient to have a user-defined array of data elements that may, or may not, be of the same data type. This is called a "Stucture" in C and a "Derived Data Type" in F90. This is a new feature that was not included in F77. Suppose, for example, we wanted to set up a data base containing information (name, age, height, weight) about 3 different people. Here is a program that uses a structure to read in this information and print the name of the lightest person: program namebase implicit none integer, parameter :: nmax=3 type person character (len=20) :: name integer :: age real :: height, weight end type person type(person), dimension(nmax) :: student integer :: i real :: weightmin do i = 1, nmax print *, 'Enter first name, age, height, weight (e.g., Bob 27 69 183)' read *, student(i)%name, & student(i)%age, & student(i)%height, & student(i)%weight enddo print *, 'Lightest student = ', student(minloc(student%weight))%name print *, 'Oldest student = ', student(maxloc(student%age))%name end program namebase The statements: type person character (len=20) :: name integer :: age real :: height, weight end type person create a defined type that can be used to name a variable later in the program. The contents are a character string (of length 20), an integer for the age, and real number for the height and weight. These contents are termed the "members" or "components" of the structure. In this example, the defined type is given the name "person" which is NOT a variable name. Rather it is a name that can later to used to define the actual variable name(s) that will have this structure. The next line creates the array student that has this structure: type(person), dimension(nmax) :: student In this line, type(person) acts just like the integer and real statements that we are used to. It can be used to define more than one variable name and these variables do not need to be arrays. Here is another example of this: type(person), dimension(nmax) :: grads, undergrads type(person) :: teacher which sets up structure arrays for grads and undergrads and a single structure for teacher. The members of the structure are referenced by appending %member to the name, where member refers to the specific member defined in the structure. So, for example, student(1)%age is the age of student(1). Note that the array index goes BEFORE the %member, not after. The % operator is also sometimes called the "component selector." Note the nifty use of the minloc and maxloc functions, avoiding the need to write do loops. This example does not really make obvious the advantages of a structure because we could easily have maintained separate arrays for name, age, height and weight. A clearer advantage of the structure approach is that we can reassign all of the members with a single statement, i.e., we can say: student(3) = student(2) rather than having to say: student(3)%name = student(2)%name student(3)%age = student(2)%age etc. However, we cannot compare structures as in the statement: if (student(3) == student(2) ) then ! ILLEGAL!! --------------------- FAST I/O in FORTRAN ----------------------- FORTRAN has generally not considered to be as flexible as C in the way that it handles i/o operations. However, with a little deviousness, it is possible to write FORTRAN programs that are capable of extremely fast i/o. The key to making FORTRAN do fast i/o is to use unformatted binary read/write operations. Suppose we have a 10000 x 10 real array that we wish to store on disk. We could output this simply as an ascii file: program testio1 implicit none real, dimension(100000,10) :: a = 1.1 integer :: i, j, count1, count2, count_rate real :: dt open (12, file='out.testio1') call system_clock(count1, count_rate) do i = 1, 100000 write (12, '(10e12.4)') (a(i,j), j=1,10) enddo call system_clock(count2, count_rate) dt = real(count2-count1)/real(count_rate) print *, 'dt = ',dt close (12) end program testio1 Notice the use of the intrinsic function, system_clock, to monitor how long the output takes. On my computer (rock, a SunBlade 100), this takes about 3.8 s and creates a 12.1 Mb file. On my Mac laptop (1.5 GHz G4) this takes 3.2 s. Next, we could output this as a binary file using similar loops over the indices: open (12, file='out.testio2', form='unformatted') do i = 1, 100000 write (12) (a(i,j), j=1,10) enddo close (12) On rock, this takes about 0.53 s (7 times faster!) and generates a 4.8 Mb file. On my Mac laptop, this takes 0.89 s (3.6 times faster). A further improvement occurs when we store the entire array directly: open (12,file='out.testio3',form='unformatted') write (12) a close (12) On rock, this generates a 4.1 Mb file in about 0.26 s. The total speed improvement compared to ASCII is about a factor of 15. On my Mac laptop, this takes 0.16 s, about 20 times faster than ASCII (The speed tests in this section were done without the -C compiler option. The -O optimizing option will produce a slight additional speed increase for testio2 and testio3) F77 yields a similar total speed improvement for testio3; however, the improvement for testio2 is not as good as in F90. Similar speed improvements can be noted upon reading these files. The message is that for speedy i/o of large data sets in FORTRAN, we should read/write large arrays in blocks. This is how Guy's GFS programs work and my routines to store travel time data. Of course, it is not always convenient to refer to data within a big array. To make it easier to keep track of things, it is often a good idea to define a data type (structure). For example, suppose our array actually consists of 10000 "lines" of data with 10 values for each line, consisting, for example, of a 4 character station name, coordinates for the station, and date/time information. In this case, we can define a suitable data type and then create an array of these data types. type station_info character (len=4) :: stname real :: slat, slon, selev, sec integer :: iyr, imon,iday, ihr, imn end type station_info type(station_info), dimension(10000) :: stlist We can now perform i/o very simply using the stlist array, while preserving the ability to access its individual parts. For example, stlist(5633)%stname is the station name of the 5633th "line" of data, stlist(232)%slat is the latitude of the 232th line, etc. F77 does not allow for structures, but a kluge is possible using the equivalence statement. For example, we could we define a small array aa with 10 elements and set these elements to be equal to the individual variable names that we want: real aa(10) common/com1/stname,slat,slon,selev,iyr, & imon,iday,ihr,imn,sec equivalence (aa(1),stname) The equivalence statement sets the values of aa to the same place in memory as the 10 variable names defined in the common block. The beauty of this is that these variables can be of mixed data types. In order to examine the record 523, we just set aa to the appropriate part of a: i=523 do 50 j=1,10 aa(j)=a(i,j) 50 continue The stname, slat, etc., are now set. This trick is used in Guy's GFS routines and in many of my I/O routines. ------------------------ ASCII vs. BINARY files ----------------------- There is a tradeoff between speed and convenience for ascii versus binary files. ascii is easier to work with and to understand and is often the preferred format when file size or speed is not a problem. However, for large data sets and data processing, binary is often better because it permits much faster I/O and the files will generally be smaller than their ascii counterparts. You are also assured of retaining the full machine precision for numerical values (you don't risk deciding later that you really should have stored 4 places to the right of the decimal, not 3). However, binary formats are somewhat less portable then ascii and require more user knowledge about their format.