The GISS GCM Good Programming Guide
Version 1.0 February 1999
Gavin Schmidt, Jeff Jonas, Jean Lerner & Ruedy Reto
The GISS GCM (and all attendant offshoots) has developed into a large
distributed effort. This guide is an effort to help integrate
GISS-wide good coding practices that improve the efficiency of the
code, make it more transparent and hopefully (in the long term) lead
to some degree of homogenization.
This guide will be split into a number of sections. Firstly, we will
highlight some of the common examples of 'bad' code and indicate some
better alternatives. These should become common sense programming
habits. The next section deals with the ways of getting rid of
unnecessary GO TO statements, and discusses more structured
approaches. We then outline some of the more useful elements of
FORTRAN 90which can be used to enhance the readability of the
code. The first appendix describes tools that are under-utilised but
can be very helpful. The other appendix contains some of the issues
that arise in porting the model code to the SGI machines. Examples are
given where we feel the correct behavior is not obvious.
This is not intended to be a comprehensive guide to FORTRAN (or to the
GISS GCM). It is intended only to highlight some of the more common
problems that occur in the model and should be seen as a useful
reference for the programmers and scientists in the building.
Please let us know about any features that could be added to this
document. They will appear in the next version (contact
gavin@giss.nasa.gov).
What to definitely avoid
- Do not calculate invariant expressions at every time step. Put
those calculations whose answer never changes in an IFIRST
section and calculate them once only at the start of every run. Note
that on the SGI machines, code segments containing an IFIRST
section must be compiled with the -static option in order to
save the local variables. However, unnecessary compilation with
-static can add a lot of overhead, slowing down the code. For
instance, QUSBxxx routines do not require any saving of local
variables and hence need not be compiled with this option. As a better
alternative, all needed local variables should be declared in a SAVE
statement and the -static option never used. Care must be
taken to make sure that all such variables have been identified before
switching over.
- Avoid calculating expressions that do not depend on the loop
variable inside the loop. Move such independent code outside the
loop.
- Do not divide! Division takes many times longer than
multiplication and so if you are dividing by the same denominator more
than once, calculate the reciprocal and multiply by that. Sometimes
the compiler can do this for you, but this cannot be relied upon.
- Organise your loops and arrays so that the innermost loop is for
the first index. For instance, in an I,J,L loop over
T(I,J,L). L should be the outermost loop, followed
by J, followed by I. Doing it any other way is
highly inefficient. If your subroutine processes things by
(I,J) grid box, and performs calculations in the vertical
(such as MSTCNV, CONDSE etc), write the internal
arrays with L as the first index.
The most egregious example of bad looping occurs right at the start of
subroutine DYNAM,
COMMON U,V,T,P,Q
COMMON/WORK6/UT(IM,JM,LM),VT(IM,JM,LM),TT(IM,JM,LM),PT(IM,JM),
* QT(IM,JM,LM)
COMMON/WORK2/UX(IM,JM,LM),VX(IM,JM,LM),TX(IM,JM,LM),PX(IM,JM)
....
DO 310 J=1,JM
DO 305 L=1,LM*3+1
DO 305 I=1,IM
UX(I,J,L)=U(I,J,L)
305 UT(I,J,L)=U(I,J,L)
DO 310 L=1,LM
DO 310 I=1,IM
310 QT(I,J,L)=Q(I,J,L)
On the face of it, there's nothing wrong. It works fine. However, if
rewritten properly the model actually runs noticeably faster. What
happened here was that the code probably started out correct just
processing UX and UT. Then someone must have added QT and decided that
they could save some coding by rearranging. (This is theory-the oldest
code that could be found was from 1985, and that has the
error. Unfortunately, MB112M9.S can't be fixed because there
are too many .U files to check for conflicts.) Anyhow, a
couple of years ago Jean went through all the dynamics routines for
the 31-layer model and rewrote all the inside-out loops. It ran 15-18%
faster when finished.
As a general rule, common blocks should not be initialised like this
since it is not transparent. See the example in the FORTRAN
90 section for the preferred method.
- Avoid exponentiation (ie. A**5). All polynomials should
be calculated using Horner's Rule (i.e. A*X**3 + B*X**2 + C*X +
D should be coded as D + X*(C + X*(B + X*A)). Factorize
as much as possible.
- Calls to mathematical functions ( log, sin, cos, etc.)
should be minimised. These are expensive. If repeated calculations are
needed, calculate them once and then save the result.
- Use constants for DO loops (not variables). Prior to the
PARAMETER statement (in FORTRAN 77), IM, JM,
LM, and related quantities were in common blocks. Now, IM,
JM, and LM are defined as parameters, and replaced by
placeholders ( IM0, JM0, LM0) in the common block. When loops
use constants, the loops are set up by the compiler, instead of on the
fly at execution time. Quantities in common blocks are treated as
variables, not constants, by the compiler so there is overhead setting
up the loop. So DO loops will run faster with IM and
JM as loop variables, rather than IM0, JM0.
However, there are still people using JMM1 and LMM1
(variables). Instead, JM-1 and LM-1 should be used
(which could not be done pre-FORTRAN 77). These are constants
and are treated as such by the compiler. As a general point, all
constants should be declared in a parameter statement.
Further information on optimization is available in the ``XL Fortran:
Optimization Guide'' handbook or the IRIX Fortran manual.
NASA GISS:
Unstructured habits Unstructured habits
Much of the GCM code (and some of the programmers!) dates back to the
punch-card days of yore. In contrast to the current fashion for 1960's
nostalgia, we believe that the code really should reflect some of the
advances made by FORTRAN 77 (and maybe even FORTRAN
90!). Many of the most confusing examples are driven by a desire
to write compact code. While this was an issue with punch-cards, it is
no-longer so important. However, if the desire for compact code is
still present, we recommend you use the new features of FORTRAN
90 which allow extremely compact code to be written in clear and
unambiguous ways (see next section).
- GOTO bad, IF (..) THEN good. Prior to the
introduction of BLOCK IF structures (and now
DO...WHILE loops) huge numbers of GOTO statements
were needed. However, there is a tendency for GOTO structures
to become very spaghetti-like, and it becomes very difficult to follow
the logic of the code. Replacing GOTO with any of the other
structures makes it much clearer for the reader (and for the
compiler).
- Initialisation across a common block using only one array is bad
habit. If the common block changes or is rearranged, this can cause
serious confusion. See the example in the FORTRAN 90 section
for a clearer way of doing this efficiently. (Also see last comment in
this section)
- Obsolete code should be deleted. Use a new version number for new
code (do not just call it 'new' or 'latest'). If you need a record,
you can always look at the previous version stored on Ipcc2 or
Ra. Commenting it out only leads to confusion. Leaving it in, and not
using it, is inefficient and confusing to the next generation of
programmers.
- Over-use of equivalence statements. While equivalence statements
are useful for writing compact code, their use should really be
restricted for operations that are applied uniformly over the
equivalenced arrays (such as initiallisation, input/output etc.)
- IF ... GOTO branches out of a DO loop. While
this is supported by our compilers and is technically correct, it is
not a recommended feature of standard fortran. Problems have occured
with compilers/optimizers that did not treat it correctly. In
particular, if a loop variable is used outside the loop subsequent to
such a branch, this sometimes gave incorrect results. This practice
should be avoided if possible. More generally, these branches
seriously inhibit proper optimization (since loop structures cannot be
changed). Use the DO..WHILE construction instead (see next
section).
- Do not use active code lines as continue statements. (This makes
Jean CRAZY!) Constructs like this
do 110 i=1,im
b = something complicated
110 a = something complicated depending on b
are not actually wrong, but they are very awkward if someone wants to
do some debugging or add an additional calculation. Then, instead of
inserting a statement, they also have to rewrite:
do 110 i=1,im
b = something complicated
a = something complicated depending on b
c = something depending on a
110 continue
This is a pain in the neck!!! It also makes it harder to see what's
going on when you do a diff on two files, because it looks like two
things were changed, not one. Please use CONTINUE or END
DO statements. (In fact with DO ... END DO, you can
dispense with line numbers entirely).
- COMMON block cautions. Two rules of clean programming
that are frequently violated in the GCM are i) out-of-bounds array
referencing (as alluded to above), and ii) declaring a common block to
be different sizes in different routines. While neither of these are
strictly illegal, they do cause problems with
optimisers/compilers. For instance, on the SGI the compiler option
-OPTreorg\_common is on by default. This pads out common
block to make them more efficient (by making sure arrays are on the
same page of memory for instance). This is now turned off in
setup. Other problems can also occur. For instance, when the
compiler encounters common blocks of different sizes it will select
the largest, except when the common block is initialised via a block
data section even if there are larger references elsewhere. If we
reduce our dependence on these kinds of violations, future problems
are likely to be minimised.
Common blocks that are used in sub-modules (such as the PBL or tracer
codes) are sometimes required in MAIN, INPUT, etc. It is much
more straightforward to put these in separate (named) files which are
then INCLUDE-d in the code. That way revisions and changes
can be quickly accommodated.
New structures from FORTRAN 90
Many of the new features in FORTRAN 90 significantly reduce
the amount of coding needed and makes the code much more
readable. There is much more to FORTRAN 90 then we can do
justice to here, but we particularly want to highlight the array
processing facilities. To use FORTRAN 90constructs, you need
only change f77 to f90 in your compilation
scripts. There is one major caveat: the GCM II' as a whole does not
run if compiled with f90. Some parts of it (the radiation?)
cannot be compiled, but others have included f90 features
with no problems. Gary's model can be run with f90. Please
ensure that any routines you modify will compile with f90
prior to including any new constructs. Please report any problems (and
solutions!) you uncover.
- Array arithmetic. In FORTRAN 90, arrays can be used
directly in an expression (as long as it is conformal) without having
to loop over the indexes. The example discussed above can be compactly
written as below with the compiler deciding the most efficient way to
loop over the variables.
COMMON U,V,T,P,Q
COMMON/WORK6/UT(IM,JM,LM),VT(IM,JM,LM),TT(IM,JM,LM),PT(IM,JM),
* QT(IM,JM,LM)
COMMON/WORK2/UX(IM,JM,LM),VX(IM,JM,LM),TX(IM,JM,LM),PX(IM,JM)
DIMENSION UALL(IM,JM,LM*3+1),UTALL(IM,JM,LM*3+1),UXALL(IM,JM,LM*3+1)
EQUIVALENCE (UALL,U),(UTALL,UT),(UXALL,UX)
....
UXALL=UALL
UTALL=UALL
QT=Q
Another example is when you wish to only set a limited number of
indexes, or ranges of the index. For this, : denotes the entire
range of an index or 2:7 for instance denotes just the range 2
until 7 (inclusive). For instance, commonly the polar boxes are set to
be the same and equal to the value at I1. This could be
written:
P(2:IM,JM) = P(1,JM)
- New constructs: CASE and DO...WHILE. The
CASE construct allows you to branch to any number of
sub-sections. It is similar to an arithmetic IF or a computed
GOTO, but is much more flexible. For instance,
SELECT CASE (ITYPE)
CASE (1) ! ITYPE 1
stuff
CASE (2:4) ! ITYPE 2,3 or 4
stuff
CASE DEFAULT ! ITYPE anything else
PRINT*,"ITYPE must be 1,2,3 or 4"
END SELECT
The DO...WHILE construct can replace convoluted constructions
involving IF and GOTO statements. In particular, it
should definitely replace constructs that include IF
statements branches out of DO loop. For example,
DO WHILE (Q(1) < QSAT)
EVAP = ....
Q(1) = Q(1) + EVAP
END DO
One last thing...
COMMENTS ARE GOOD! Please, please, please comment any code you
write. It doesn't need to be very verbose, but important calculations,
break points, branches etc. should be flagged. At some stage, other
people will try to work through it - please think of them! In
particular, messages printed out before a stop due to unrealistic
conditions should be verbose rather than cryptic.
Appendix I: Helpful tools
The dbx debugger
dbx is a powerful debugging utility that can be used quickly
to diagnose some common problems. If you have a core dump, then go to
the run directory and type
dbx program.exe core
At the prompt, typing 'where' will give the routine and line number of
where it died. Compiling with the -g option increases the
amount of information available, i.e. 'print variable' will display
the value of the variable at the time of the core dump.
The KAP preprocessor
A Preprocessor rewrites your fortran code before sending it to the
compiler and can be very useful in identifying bottlenecks and
inefficiencies in your code. Using the preprocessor can lead to
significant improvements in time (average 10-15% faster for Gary's
model). It is called via an option to the fortran compiler but is not
available on all machines.
Appendix II: Modifications necessary for SGI computers
- All FORTRAN and C utilities have to be recompiled.
- Most ksh scripts are ok, some have to be modified - test
them. Please be sure to add #!/bin/ksh at the top of all
scripts, so that those who have different login shells can use them.
- xlf .... has to be replaced by f77 .... (or f90) There has to be a space between -o and
name-of-executable-file
OPTIONS:
- in model routine compilations the options -n32 -static or
-64 -mips4 -static are needed (the 2nd gives better
performance).
- -static guarantees that local variables in a subroutine
are saved after returns, i.e. IFIRST-constructions should work.
- -old_rl should be used with direct-access files to make it
compatible with the IBMs
- -OPT:reorg_commonoff is needed with the new release
- -mp (in the link step) will allow large jobs
to run much faster, if your .profile contains the line:
export PAGESIZE_DATA64; export PAGESIZE_STACK64; export PAGESIZE_TEXT64
- -Wl,-woff,134 -Wl,-woff,15 arguements for f77/ f90 used when linking (in setup or mkexe)
removes the pointless warnings, but leaves any important ones.
- Optimizations
( -O3 cannot be used for individual subroutines for -n32 )
what worked so far is:
f77 -n32 -static -O2 -c xyz.f (-> xyz.o) or
f77 -64 -static -O3 -c xyz.f (-> xyz.o)
for more speed:
f77 -64 -mips4 -static -O3 -OPT:fold_arith_limit1409
-O2 is almost as fast and safer than -O3 (optimization level)
link options:
-mips4 lfastm (uses library of fast math.routines)
- Differences in FORTRAN (missing/obsolete features etc):
- READ(...,NUMnbytes) is unavailable
- BLOCK DATA have to be named
- 'open' cannot be the name of a COMMON BLOCK
- T,F cannot be used to initialize logical variables, unless
PARAMETER (T.TRUE.,F.FALSE.) is added
- Use STOP rather than RETURN in MAIN.
- The function ERFC only takes real*4 args, use DERFC
for real*8 args.
- The compiler warns if you use CALL SUBR(A,..) where A is a
variable whereas A is an array in the subroutine (using of
course only A(1))
- If the arrays and their dimensions are passed to a subroutine,
the dimensions have to be passed also in each entry that
passes the arrays (the original HNTRP does not work on
the SGIs)
- The following construction does NOT work (with or w/o -static)
SUBROUTINE ABC (v1,...)
(no v1... statement)
ENTRY XYZ (v1 is NOT an argument)
expression involving v1 (result unpredictable if
reached after CALL XYZ)
Quick (but ugly) fix:
SUBROUTINE ABC (v1xxx,...)
v1v1xxx
- NAMELIST READ('string',NML...) is unavailable ; namelists have to be read
from external files: back to 66-version (using unit 8 instead of 14)
- System calls to save rsf, acc files did not work reliably. Use
simple OPEN instead as in model II' ( MB112M9.S).
- Model-specific changes:
- remove all READ(....NUM...) instances ( INPUT,readt
in UTIL..)
- replace internal NAMELIST reads by reads from unit 8 ( INPUT,2x)
- name all BLOCK DATA (diagnostics, rad.routines, input,soils)
- replace RANDU-subroutine (add RANDSGI)
