An Introduction to Fortran

1
An Introduction to Fortran

• Lecture 2
• Mar 20 2009
• Kevin Keay

2
Outline

• Lecture 1 10-12 PM (Thu)
• Lecture 2 10-12 AM (Fri)
• Lab session 1 2-330 PM (Thu)
• Lab session 2 At your leisure

3
Arrays

• See also Primer (Mardling 1997)
• A one dimensional array x is the representation
  of a vector x in mathematics
• A vector x (x1 x2 xn)T
• is represented as
• real x(100)
• where we are assuming n 100
• or
• integer maxn
• parameter (maxn100)
• real x(maxn)
• We can process x for subscripts 1 n where n
  maxn

4

• X1 ?? x(1)
• X2 ?? x(2)
• Xn ?? x(n)
• The individual elements of the array x may be
  treated like simple variables
• (1) y x(2) c
• (2) do k3,n-1
• z(k) x(k) r 2.518
• enddo

5

• An implied DO loop is useful for reading arrays
• read(1,(2X,F9.4))(y(k),k1,n)
• (we assume that n is known)
• Arrays can be of any type
• (1) integer count(50)
• (2) integer maxm
• parameter (maxm1000)
• character20 name(maxm)
• i.e. name(k) is a text variable of up to 20
  characters for k1 maxm

6
2D arrays

• A 2D array is the representation of a matrix in
  mathematics. For a r rows x c columns matrix x
• xij ?? x(i,j) i1 r, j1 c
• integer maxr,maxc
• parameter (maxr50,maxc10)
• real x(maxr,maxc)
• The maximum number of rows is maxr
• The maximum number of columns is maxc
• See Primer (Mardling 1997) Section 3.2

7

• A pair of nested DO loops is useful for
  processing 2D array elements
• do j1,nlat
• do i1,nlon
• y(i,j) z(i,j) a b
• enddo
• enddo
• Here i corresponds to longitude and j to
  latitude
• Note the array element z(i,j) may be treated as
  a simple variable

8

• Evaluate y(i,j) z(i,j) a b
• This is effectively what occurs
• j1 is set by do j1,nlat
• do i1,nlon
• y(i,1) z(i,1) a b
• enddo
• j2 is the next value of j set by do j1,nlat
• do i1,nlon
• y(i,2) z(i,2) a b
• enddo
• jnlat is the last value of j set by do j1,nlat
• do i1,nlon
• y(i,nlat) z(i,nlat) a b
• enddo

9

• An 2D implied DO loop is useful for reading and
  writing 2D arrays
• read(1,)((y(i,j),i1,nr),j1,nc)
• For nr2 and nc3 this is equivalent to
• read(1,)y(1,1),y(2,1),y(1,2),y(2,2),
• y(1,3),y(2,3)
• Note Work from the outside in
• i.e. j is set first, then i goes from 1 to nr
• Furthermore, the 2D array is read in column
  order
• i.e. by j.

10
Binary (unformatted) files

• The files encountered in the primer and reference
  are formatted files that are appropriate for text
  data. We open such files with a statement like
• open (1,filepress.dat)
• read (1,)x
• Much of the output from GCMs and related
  reanalysis projects e.g. NCEP is in a binary form
  i.e. not readable to the naked eye, like text.
  Such files are opened and read in a different
  way. Assume we have a binary version of press.dat
  called pressb.dat
• open (1,filepressb.dat,formatunformatte
  d)
• read (1)x
• In the open statement we need the extra keyword
  format set to unformatted that means binary in
  Fortran. Also in the read statement we just have
  the unit number without any format i.e. (1)

11
Conmap (CMP) format

• We have software that takes the binary files from
  (say) NCEP Reanalysis and converts them to a
  simple binary format which we call conmap (CMP)
  format. This is capable of representing a
  variable e.g. pressure, on a longitude-latitude
  grid.
• A fragment of code to read a CMP file) is

12

• implicit none
• real xmiss
• parameter (xmiss 99999.9) ! conmap
  missing value code
• integer num
• parameter(num 200) ! Max. size of
  lat.-lon.
• ! grid
• real xlats(num),xlons(num) ! Arrays for
  lats and lons
• real dat(num,num) ! Array for
  input data
• character80 head1 ! Header
  (description)
• character80 infile ! CMP file name
• open(1,fileinfile,form'unformatted')
• read(1)nlats
• read(1)(xlats(i),i1,nlats)
• read(1)nlons
• read(1)(xlons(i),i1,nlons)
• read(1)head1
• read(1)((dat(i,j),i1,nlons),j1,nlats)
  !See conmap.f
• close(1)
• write(,'('' No. lats, no. lons
  '',2I6)')nlats,nlons

13

• We use a parameter statement to set the maximum
  size of our longitude and latitude arrays to 200
  points (num). Thus our two-dimensional data array
  (matrix)
• named dat is a grid of maximum size 200 x 200
  points
• The integers nlons and nlats hold the actual
  number of longitudes and latitudes e.g. nlons
  144, nlats 73
• The array xlons holds the actual longitude values
  e.g. 0.0, 2.5, and the array xlats holds the
  actual latitude values from
• south to north e.g. 90.0, -87.5,
• The array xlats is read using an implied DO loop
• read(1)(xlats(i),i1,nlats) is
  equivalent to
• read(1)xlats(1),xlats(2),,xlats(nlats)
• The array xlons is read in a similar way. There
  is an 80 character text header (head1) containing
  some information about the data.

14

• Finally we read in the two-dimensional array dat
  using a pair of nested implied DO loops
• read(1)((dat(i,j),i1,nlons),j1,nlats)
• ? i is longitude
  index and j is latitude index
• This works from the outside to the inside
  loop. We start with j 1 and read the inner loop
  (dat(i,1),i1,nlons). Then j 2 and we read
  (dat(i,2),i1,nlons) and so on.
• For each latitude index j we read the values of
  dat at each longitude index i i.e. we read the
  data on latitude circles from south to north. If
  we specified each array element individually we
  would need
• read(1)dat(1,1),dat(2,1),,dat(nlons,1),
• dat(1,2),dat(2,2),
• dat(nlons,2),,dat(1,nlats),dat(2,nlats),dat(n
  lons,nlats)

15

• To write out a CMP file we just replace read by
  write in the above code.
• The CMP format has a code to represent missing or
  undefined data (99999.9) as set by xmiss. This is
  useful to process files with missing data since
  we want to exclude these from our summations etc.
  See statcon.f .
• If we expect missing values then we should check
  for them.
• This piece of code computes the average of
  non-missing (live) values at each gridpoint
• do j1,nlat
• do i1,nlon
• if (dat(i,j).eq.xmiss)then
• nmiss nmiss 1 ! Count of missing
  values (for information only)
• else
• nsum(i,j) nsum(i,j) 1 ! Count of live
  values at this gridpoint
• sum(i,j) sum(i,j) dat(i,j) ! Sum of
  live values
• endif
• enddo
• enddo
• do j1,nlat

16
Command-line arguments getarg, iargc

• These are useful intrinsic (built-in) routines to
  simplify
• input and output.
• Imagine we have a program called jason1. By using
  the above routines we can read filenames and
  other arguments from the command line. In our
  case we want to give
• an input file (argument 1) and an output file
  (argument 2)
• jason1 press.dat mean.dat
• ?argument 1 ?argument 2
• A program fragment to do this is

17

• character optarg80, infile80, outfile80
• integer narg
• narg iargc() ? This will be 2 in the above
  example i.e. 2 arguments
• write(,)No. of arguments ,narg
• call getarg (1,optarg) ? Get first argument
  from command line and put into optarg
• infile optarg ? This variable holds the name
  of the input file
• call getarg (2,optarg) ? Get second argument
  from command line and put into optarg
• outfile optarg ? This variable holds the name
  of the output file
• open(1,fileinfile) ? Open the input file
• open(2,fileoutfile) ? Open the output file

18

• We can also read numerical arguments from the
  command line
• jason1 press.dat mean.dat 24.5
• ?argument 1 ?argument 2
  ?argument 3
• Argument 3 can be read by adding
• real x
• call getarg (3,optarg) ? Get third argument from
• 
  command line and put into optarg
• read(optarg,)x ? Do an internal read from
  optarg i.e. read x
• from optarg
  
• write(,)x ,x

19

• This approach is very useful for handling a large
  set of input files
• jason1 press.9601?? mean.dat
• Under UNIX or Linux the argument press.9601?? is
  expanded into a set of filenames which differ in
  the last two characters. Note that at least one
  of these files must exist. Assume that the files
  are called press.960101, press.960102,
  press.960103, , press.960107 i.e. daily pressure
  files for Jan 1 1996 to Jan 7 1996. Then the
  command line will be expanded as if the
  individual files were given i.e. as if we had
  specified
• jason1 press.960101 press.960102
  press.960107 mean.dat
• ?Argument 1
  ?Argument 2 ?Argument 7
  ?Argument 8
• Such an approach is used in statcon.f (see Lab
  Session
• notes example 9).

20
Example 1 Calculation of the mean of a set of
numbers

• Open a text file containing a single column of
  numbers ( 1000)
• Read in the numbers (1 per line) into a real
  array x this process also gives the sample size
  n
• Compute the mean meanx of the array x
• Write the mean and sample size on the screen
• Note the use of the continuation character a
  in column 6

21

• program mean
• c
• c Purpose Reads in a textfile
• c containing one value per line (x),
• c computes the mean of the data,
• c and writes the mean value on the
  screen.
• c
• c Author Kevin Keay Date 3/3/03
• c
• implicit none
• c
• character80 infile
• integer i,n
• real x(1000) ! Array of size 1000
• real meanx
• c
• write(,)'Enter input file name'
• read(,'(A)')infile
• c

22
Example 2 Calculation of the mean of a set of
numbers using a subroutine and a user function

• The code to compute the mean is relegated to a
  subroutine named getmean
• Alternatively, a user function getmean2 may be
  used

23

• program mean2
• c
• c Purpose Reads in a textfile
• c containing one value per line (x),
• c computes the mean of the data,
• c and writes the mean value on the
  screen.
• c Uses a subroutine meanx to compute
  mean.
• c Author Kevin Keay Date 3/3/03
• c
• implicit none
• c
• character80 infile
• integer i,n
• real x(1000) ! Array of size 1000
• real meanx, meanx2
• c
• c User functions
• c
• real getmean2

24

• c-------------------------------
• subroutine getmean (n,y,meany)
• c
• implicit none
• integer n
• real y(n)
• real meany
• integer j
• c
• meany 0.
• do j1,n
• meany meany y(j)
• enddo
• meany meany/n
• write(,)'Subroutine getmean - meany
  ',meany
• c
• return
• end
• c-------------------------------

25
Important points

• The subroutine name (getmean) does not have a
  type but the user function (getmean2) does (real
  function)
• Note that names of routines in a program must be
  unique
• The parameters in the main program (mean2) and
  the subroutine (getmean) must match
• i.e. have a one-one correspondence
• mean2 call getmean (n,x,mean)
• getmean subroutine getmean (n,y,meany)
• n ?? n (integer)
• x ?? y (real array of size 1000)
• mean ?? meany (real)
• In the main program real x(1000)
• In the subroutine real y(n)
• Technically we should have y(1000) but y(n) is
  used with the implication that n1000

26

• The user function returns the mean via its name
  (getmean2) i.e. the function is a special kind of
  variable
• Like a subroutine the parameters must match
• In this case the same names are used in the main
  program and the function
• As in the subroutine we could have used
  different names but they have to match in terms
  of type

27
Example 3 conmap (CMP) format

• This illustrates how to read a CMP file (binary)
  and write it out in a readable format

28

• program readcon
• c
• c Purpose Reads in a conmap file and outputs
  in ASCII (readable) format
• c to a text file.
• c Write out in two ways - compact and
  verbose.
• c
• c Author Kevin Keay Date 3/2/99
• c
• implicit none
• c
• real xmiss
• parameter (xmiss 99999.9) ! conmap
  missing value code
• integer num
• parameter(num 200) ! Max. size of
  lat.-lon. grid
• real xlats(num),xlons(num) ! Arrays for
  lats and lons
• real dat(num,num) ! Array for
  input data
• c
• character80 infile ! Name of
  input files
• character60 outfile ! Name of
  output file

29

• c
• c Read in conmap data
• c
• open(1,fileinfile,form'unformatted')
• read(1)nlats
• read(1)(xlats(i),i1,nlats)
• read(1)nlons
• read(1)(xlons(i),i1,nlons)
• read(1)head
• read(1)((dat(i,j),i1,nlons),j1,nlats) !
  See conmap.f
• close(1)
• c
• c Write header to screen
• c
• write(,'('' No. lats, no. lons
  '',2I6)')nlats,nlons
• write(,'(1X,A)')head
• c
• c Write out conmap data in ASCII format
• c Do it two ways - compact and verbose

30
Important notes

• For this binary file we use read(1) not read(1,)
  or a format
• The parameter statement is used to set the values
  of xmiss and num (if num needs to be larger then
  only one change is required).
• The intrinsic (built-in) function iargc() is used
  to detect how many arguments (items) are typed on
  the command line after the program name (readcon)
• If there are no arguments i.e. just readcon
• then a usage message is written and the
  program is forced to finish with stop
• If there are not exactly two arguments then
  this is an error and the program is forced to
  finish too

31

• Otherwise we use the intrinsic subroutine getarg
  to read each of the two arguments that were input
  from the keyboard (the first is the input CMP
  file, the second is the name of the output text
  file)
• We then open the CMP file (unit 1) and read in
  all of the map information i.e. header, lat-lon
  values and the gridded data
• The map information is the written to the output
  text file (unit 2) in two ways (compact and
  verbose)

32
Extensions

• We can manipulate the gridded data in the
  two-dimensional array (matrix) named dat
• Note that i is longitude index and j is
  latitude index
• xlons(i) is the longitude value and xlats(j)
  is the latitude value
• of the gridpoint dat(i,j)
• For instance the data could be temperature in C
  and we want Kelvin
• do j1,nlats
• do i1,nlons
• dat(i,j) dat(i,j) 273.15
• enddo
• enddo
• In this case dat is being overwritten but we
  could
• introduce an extra array e.g. real tempk
  (num,num)
• and use
• tempk(i,j) dat(i,j) 273.15

33

• We can output this new gridded data to an output
  CMP file
• character80 ocmpfile
• c
• c Write out conmap data
• c
• open(2,fileocmpfile,form'unformatted')
• write(2)nlats
• write(2)(xlats(i),i1,nlats)
• write(2)nlons
• write(2)(xlons(i),i1,nlons)
• write(2)head
• write(2)((tempk(i,j),i1,nlons),j1,nlats)
• close(2)
• We could have a new header as well e.g. ohead,
  and set it to reflect the change to Kelvin
• We can reuse unit 2 as long as it is not already
  open

34
Example 4 Compute the correlation of two series

• This illustrates multiple subroutines and a old
  user function (probst) from a journal (1968)
  (this has capitals!)
• The program has a nice modular design

35

• program corr2
• c
• c Purpose Computes the correlation
  coefficient of two series x and y of length n.
• c The (x,y) pairs are given in a text file.
  It is assumed that there are no
• c missing values.
• c The correlation and its significance are
  printed on the screen.
• c Note The statistical test is H0 corr 0
  v H1 corr ltgt 0
• c i.e. a two-sided test.
• c As part of the computations the means and
  standard deviations of the two series
• c are also produced.
• c
• c Author Kevin Keay Date 3/2/99
• c
• implicit none
• c
• integer maxn
• parameter (maxn1000) ! Max. no. of (x,y)
  pairs
• real x(maxn),y(maxn)
• c

36

• c
• c ---------------------------------------
• c
• subroutine readdata (infile,n,x,y)
• c
• c Purpose Reads in (x,y) pairs from a
  textfile.
• c
• c Author Kevin Keay date 3/2/99
• c
• implicit none
• c
• integer maxn
• parameter (maxn1000) ! Max. no. of (x,y)
  pairs
• real x(maxn),y(maxn)
• c
• character80 infile
• integer i,n
• c
• c Read (x,y) pairs

37

• c ---------------------------------------
• c
• subroutine corr (infile, n,x,y,corrxy,pval)
• c
• c Purpose Computes the correlation between x
  and y.
• c Also returns the p-value (significance) of
  the correlation.
• c
• c Author Kevin Keay date 3/2/99
• c
• implicit none
• c
• integer maxn
• parameter (maxn1000) ! Max. no. of (x,y)
  pairs
• real x(maxn),y(maxn)
• c
• character80 infile
• integer n
• real corrxy,pval

38

• c
• c ---------------------------------------
• c
• integer function ilen(string)
• c
• c Function ilen
• c Purpose Returns the position of the last
  non-blank character
• c in string. If string is entirely
  blank ilen -1 .
• c
• implicit none
• character() string
• integer i
• do ilen(string),1,-1
• if(string(ii).ne.' ')then
• ilen i
• return
• endif
• end do
• ilen -1 ! string is entirely blank

39

• REAL A, B, C, F, G1, S, FK, T, ZERO, ONE,
  TWO, HALF
• C
• C G1 IS RECIPROCAL OF PI
• C
• DATA ZERO, ONE, TWO, HALF, G1
• /0.0, 1.0, 2.0, 0.5, 0.3183098862/
• C
• IFAULT 1
• PROBST ZERO
• IF (IDF .LT. 1) RETURN
• C
• IFAULT 0
• F IDF
• A T / SQRT(F)
• B F / (F T 2)
• IM2 IDF - 2
• IOE MOD(IDF, 2)
• S ONE
• C ONE

40

• The examples in Example 11 of the lab session
  illustrate several ways to put together this
  program
• The data used are the two time series shown in
  the following graph

41
NCEP Reanalysis extratropical cyclone count for
30-50S and SH temperature
42
(No Transcript)