Numerical Weather Prediction (NWP), for both operational and research purposes, requires only fast computational speed but also large memory. A technique for solving the Primitive Equations for atmospheric motion on the CYBER 205, as implemented in the Mesoscale Atmospheric Simulation System, which is fully vectorized and requires substantially less memory than other techniques such as the Leapfrog or Adams-Bashforth Schemes is discussed. The technique presented uses the Euler-Backard time marching scheme. Also discussed are several techniques for reducing computational time of the model by replacing slow intrinsic routines by faster algorithms which use only hardware vector instructions

Tuccillo, J. J.

English

NASA Technical Reports Server

Memory Efficient Solution of the Primitive Equations 
for Numerical Weather Prediction on the CYBER 205 
James J. Tuccillo 
Systems and Applied Sciences Corporation 
5809 Annapolis Road 
Hyattsville, MD 20784 
1. INTRODUCTION 
Numerical Weather Prediction (NWP), for both operational and research purposes, 
requires not only fast computational speed but also large memory. In this paper I will 
discuss a technique for solving the Primitive Equations for atmospheric motion on the 
CYBER 205, as implemented in the Mesoscale Atmopsheric Simulation System (MASS) 
(Kaplan et. al., 19821, which is fully vectorized and requires substsntially less memory 
than other techniques such as the Leapfrog or Adams-Bashforth Schemes. The technique 
to be presented uses the Euler-Backard time marching scheme. 
Also to be discussed will be several techniques for reducing the CPU time of the 
model by replacing %lowt’ intrinsic routines by faster algorithms which use only hardware 
vector instructions. 
35 
https://ntrs.nasa.gov/search.jsp?R=19840012145 2020-03-22T08:32:32+00:00Z
2. MODEL BACKGROUND 
2.1 Description 
MASS is a hydrostatic primative equation model which is run over a limited 
area. The model forecast the 3-dimentional structure of wind, pressure, 
temperature and moisture. The actual domain of coverage, along with the 
horizontal distribution of grid points, is depicted in Fig. 1. The characteristics of 
the model are listed in Table 1. 
2.2 Uses and Support 
The model has been applied primarily to the problem of forecasting the 
atmospheric environment within which severe local storms (severe thunderstorms 
and tornadoes) are likely to develop. It has also been applied to the problems of 
forecasting and investigating east coast cyclogenesis, upper level turbulence and 
shear, and boundary layer transport. Support for the model development has been 
provided by NASA/Goddard using the computational facilities of NASA/Langley 
(CYBER 203) and NASA/Goddard (CYBER 205) 
2.3 History 
The original version was implemented on a 500K word CDC STAR 100 Vector 
Processor at NASA/Langley in the late 70’s using 64-bit FORTRAN. The 
availability of the SL/l programming language at Langley, which permitted easy 
access to the 32-bit instruction set on the STAR 100, resulted in an effective 
doubling of the memory and the model was recoded with larger vectors. This 
allowed for an increase in the area over which the model was run while maintaining 
the same horizontal and vertical resolution. 
36 
1 
5 
9 
13
 1
7 
21
 2
5 
29
 3
3 
37
 4
1 
‘IS
 4
9 
53
 5
7 
61
 6
5 
69
 X
I 
77
 6
1 
85
 6
9 
93
 9
71
01
10
51
09
11
31
17
12
1 12
5 
83
 
89
 
K 61
 
77
 
n 88
 
8s
 
61
 
!i7
 
53
 
¶S
 
45
 
41
 
37
 
93
 
28
 
26
 
21
 
17
 
13
 
8 5 1 
1;
 
1.
1’
 
=.
. 
. .
 . 
. .
 . 
. 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
-..
.;.
 
.b
 . 
. .
 . 
. 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
. .
 . 
. .
 
J 
5 
8 
13
 1
7 
21
 2
6 
29
 3
3 
37
 4
1 
45
 *
II 
53
 6
7 
61
 6
5 
69
 7
3 
-K
’ 8
1 
85
 8
9 
93
 8
71
01
10
61
09
11
31
17
12
11
25
 
FI
G
, 
1 
D
O
M
AI
N
 
O
F 
C
O
VE
R
AG
E 
BY
 
M
AS
S M
O
D
EL
 
AN
D
 T
H
E 
H
O
R
IZ
O
N
TA
L 
O
F 
G
R
ID
 
PO
IN
TS
 
D
IS
TR
IB
U
TI
O
N
 
TABLE 1 CHARACTERISTICS OF MASS MODEL 
0 
0 
0 
0 
0 
0- 
0 
0 
0 
0 
m(DEsCF?lPTIoN) 
Hvmmnc fhHITM hMNI& 
TERAIN FQLOWINGSIGHA-P Ctmrmm 
LIMITED Am DOHAIN 
CM~IMGRIDU~A P~SE~~PHIC~~AP(ARAKAYA"A" Grub> 
~TH OFWR Accmm Hmzmm PACE DFFE~NCING 
&OmmAccumVm~~ SPACE DFFGNING 
~NDOIW Accum TIHE DIF~~~NG 
MIADATAI~DERIVED~~H THE LFtl Arwx ~RAWINSO~~DES 
INmnrZA~acsBAsEDOIJMECACULGaCV~ATI~ 
PHvsrcS 
-LMGESCMf%fCIPIlATION 
-~ETMYB'QWJD~YLAY~ 
- DRY CoNwpIbN 
* - Mosr C~NVECTI~~J(UPIDERDEV~OPMENT) 
50 Kn GRIDSPACINGAT 45O11 
19Ewuv SPACEDLAY= 
128 X 96 COHPUTATIONL DOWN 
TI HE DEPENDENT MUNDARY Co~~rnws 
COMFWEHEHSIVE~NTEMCTIVE DIAGWJIC PACKA~CE ON THE FRONT END 
-vEKTIc#PRoms 
- vEmc# CRassEcTIoNs 
-Cmsrr\m P~ssla SWACES 
- TIHE Hmm~ 
- TFWXXIES 
- VFWKATIONSTAIKTICS 
38 
In the spring of 1980, the STAR 100 was upgraded to a lm word CDC CYBER 
203. The new machine effectively had twice the memory of the STAR 100. The 
area over which the model is run was again expanded and the vertical resolution 
was increased from 12 to 14 vertical layers. 
In the spring of 1983, the model was transferred to the NASA/Goddard 
CYBER 205. The model was recoded in CDC FORTRAN 2.0 using 32-bit 
arithmetic. After being successfully benchmarked against the Langley version, the 
vertical resolution was again increased from 14 to 19 layers. The Goddard version 
of MASS on the CYBER 205 executes approximately 3 times faster than the 
Langley version on the CYBER 203. This can be explained by 
1)Reduction in cycle time from 40 to 20 NS. 
2)Linked triad instruction on the CYBER 205. 
3)Faster gather/scatter instruction. 
4)Coding differences. 
3. EQUATION SET 
The model utilizes a standard primitive equation set cast in a terrain followingrp 
coordinate system. As indicated earlier, the forecasted variables are the 3-D 
distribution of wind, pressure, temperature and moisture. The basic prognostic equations 
are given below where u and v are x and y coordinate momentum, T is temperature, q is 
the moisture mixing ratio andnis the pressure at the terrain minus the pressure at the 
top of the model. 
39 
Three diagnostic equations close the system and are given below where & is the 
vertical velocity, 4 is the geopotential energy andCc)is the vertical velocity in pressure 
coordinates. 
The boundary conditions are 
and the definitions for and Vare 
o-= P-Ptop 7-l-= P SJR - P t OP 
the remaining variables are 
40 
m= mapscale grid transformation factor 
c= P specific heat at constant pressure 
R= gas constant for dry air 
P sup= pressure at the terrain 
Ptop= pressure at the top of the model 
X = horizontal eddy diffusivity 
4. GRID SYSTEM 
The technique for solving the differential equations is to discretize the equations 
into finite difference form and solve them on a 3-D grid. The horizontal grid employed is 
the Arakawa “A” grid where all dependent variables are defined at all grid points. The 
vertical grid is staggered so that u, v, T and q represent layer averages defined at the 
midpoint of each layer and and are held at the layer interfaces. The third diagnostic 
variable, w, is held with u, v, T and q. This structure is represented in Fig. 2. 
5. NUMERICAL TECHNIQUE 
5.1 Horizontal Space Derivatives 
The fourth order accurate finite difference approximation to an x-direction 
space derivative for an arbitrary variable y/ is given below 
where i is a horizontal index; An analogous formula is used for y - direction 
derivatives. 
41 
s-- lJ,V,T,Q,w - - - K+l 
--- U,V,TJ,w - - - K 
64 K-X 
--- U,V,T,Q,o - - - K-l 
l 
FIG, 2 VERTICAL GRID SYSTEM OF MASS 
42 
- 
5.2. Vertical Space Derivatives 
A second order accurate finite difference formula is used to approximate the 
vertical advection terms of the u,v, T and q prognostic equations. The 
representation, for an arbitrary variable p, is given below 
where k is a vertical index. 
5.3 Time Derivatives 
A second order accurate approximation to the time derivatives is used. The 
Euler-Backward Technique has the properties of frequency dependent damping and 
no computational mode. For an arbitrary variable the finite difference 
representation is given as 
Prediction 
Correction 
where n is a time level index and * refers to a intermediate time level. 
43 
This scheme requires the storage of only one time level of information (time 
level n) whereas other explicit schemes such as the Leapfrog Scheme requires the 
storage of at least two time levels (n and n-l). The penalty is that twice the 
computational work is required as compared with the Leapfrog scheme. 
6. BASIC MEMORY REQUIREMENTS 
As mentioned earlier, the Euler-Backward scheme for time marching the 
prognostic equations for the 3-D structure of wind, pressure, temperature and moisture 
requires the storage of only one time level of information. The * ‘ed time level is an 
intermediate time level and only needs to be as deep (with respect to the vertical) as is 
required to solve the equations at a layer. It should be noted that only the vertical 
advection terms couple the model layers together and that to solve the equations at leyer 
k requires the dependent variables at layers k+l, k and k-I. Therefore, the * ‘ed time 
level only needs to be 3 deep (it holds the prediction values to be used during the 
correction step) and can be reused for the solution of each layer. 
Given that the 19 model layers contain 128 x 96 grid points each, the basic 
memory required is 
u (128, 96, 19) 
v (128, 96, 19) 
T (128, 96, 19) 
q (128, 96, 19) 
pi (128, 96) 
ustar (128, 96, 3) 
vstar (128, 96, 3) 
44 
tstar (128, 96, 3) 
qstar (128, 96, 3) 
pistar (128, 96) 
If an additional layer were to be added only the u, v, T and q arrays would be 
increased. The ustar, vstar, tstar and qstar arrays are always dimensioned 3 deep and 
this is a function of the vertical advection terms which require 3 layers of storage to 
solve the equations. 
In contrast, the Leapfrog scheme would require 2 sets of arrays dimensioned 128 x 
96 x 19, therefore, there is a considerable memory savings with the Euler-Backward 
Scheme. A technique developed by Tuccillo (1983) shows some promise in reducing the 
computational work by increasing the premissable timestep. 
7. METHOD OF SOLUTION 
The method of solution is depicted in Fig. 3 and shows the sequence of steps 
required to solve the equations at all layers. Prediction is the step that advances the 
solution from the n to the * time level and correction is the step that advances the 
solution from the * to the n+l time level. It there are NZ layers then there are 2*NZ 
number of steps required to advance the solution one time step. The number above each 
line represents the order of solution where the first step is to perform prediction for 
layer 1, the second step is prediction at layer 2, the third step is correction at layer 1 
and so on. After correction (the 2*NZ step) at layer NZ is finished the solution has been 
advanced one time step. 
45 
LAYER 
NZ 
2 + NZ-2 2*NZ 
NZ- 1 
2 + NZ-4 2*NZ-1 
3 
2 
1 
2 l NZ-3 
l 
l 0 
0 
4 
0 
2 5 
1 3 
PREDICTION CORRECTION 
FIG, 3 SEQUENCE OF STEPS TO ADVANCE THE SOLUTION 
ONE TIMESTEP 
46 
The *Ied arrays are reused for each layer and the calculations for each layer are 
fully vectorized where the vector lengths are NX*NY or 12288. For this vector length 
the machine is computing at about 98% of its maximum rate. 
8. BOUNDARY CONDITIONS 
Since MASS is a limited area model, as opposed to a global model, the solution at 
the horizontal boundaries needs to be specified. The technique for specifying the 
boundary conditions consist of blending externally calculated values using a weighted 
average formula which is represented by 
where W = 0 on outer column and row 
w= 0.333 on first column and row in 
W = 0.666 on second column and row in 
w= 1.0 on third column and row in 
It should be pointed out that this technique produces an overspecification at the 
boundary and higher horizontal diffusion is required near the boundaries to control noise 
generation. 
This technique is vectorized by holding the externally specified boundary 
tendencies in a vector and using the scatter instruction to expand them into the correct 
positions prior to computing the weighted average. This technique minimizes to amount 
of storage required. 
47 
9. PROGRAMMING TECHNIQUES 
The code Is completely vectorized in the horizontal. The average vector length 
is about 12000 which represents the number of horizontal grid points. There is a loop 
over the vertical layers. 
Some specific techniques used during the coding are 
0 32-bit arithmetic 
Sensitivity tests have indicated that 32-bits provides enough precision. 
Using 32-bits effectively doubles the real memory and halves the execution 
time. 
0 Explicitly Vectorized 
The code does not depend on automatic vectorization by the compiler. 
All descriptors are set up with DATA and ASSIGN statements. Special Q8 
calls are used where required. 
.o Diadic and Triatic Structure 
All vector statements are written in a diadic structure (triadic when 
linked triads are created) to minimize compiler generated dynamic space 
which may cause paging. 
0 Subroutines are kept small enough so that the Register File is not 
overflowed. 
Subroutines which have more local variables then the size of the 
register file (approximately 200) can be inefficient since loads from 
memory must be executed. AR subroutines are kept small enough so that 
the swap instruction can load all necessary local variables at entry. 
48 
0 Parameter Statement Used for V.ector Dimensions 
Vector dimensions are easily changed by changing parameter values.. 
0 Factoring of Equations to yield Linked Triads 
The sequence of instructions have been arranged to yield the maximum 
number of linked triads. 
0 Run Only in Real Memory 
No page faults are generated during the interatlve time marching. 
0 Vectors are Grouped on Large Pages 
All large vectors are placed in common and grouped on large pages 
using loader options. 
0 Bit Vectors vs. Gather/Scatter 
For those situations where control store or gather/scatter can be 
applied, an analysis using the nominal performance figures for each 
instruction was performed and the most CPU or memory efficient 
techniques was applied. 
10. TECHNIQUES FOR REDUCING CPU TIME 
A 24-hour simulation with the model requires 1312 timesteps. Each timestep 
requires the evaluation of 2*NZ natural logs (for 12288 grid points). This required 
approximately 22 mins of CPU time using the 32-bit FORTRAN VHALOG function. 
Since the range of arguments for the natural log function was known, a more efficient 
49 
technique was incorporated where the natural log was approximated with a series 
factored using Homer’s Rule. The evaluation requires 11 vector instructions, nine of 
which are linked triads, and runs approximately 40 times faster than the FORTRAN 
intrinsic function. This technique reduced the CPU time spent evaluating natural logs to 
30 sets. 
Other techniques for reducing CPU time consist of approximating the ** 
FORTRAN function with series of square roots (square root in a hardware instruction) 
and inverting scalars to generate vector multiplies instead of vector divides. 
The version of MASS implemented on the CYBER 205 at NASA/Goddard requires 
13 large pages of memory and 15 minutes of CPU time (same as wall time) for a 24 hour 
simulation over the area depicted in Fig. 1. 
il. EXAMPLE OF OUTPUT 
MASS at Goddard features a comprehensive postprocessing system to produce 
output from the model for interpretation. The post processing system runs interactively 
and produces hard copies on a GOULD electrostatic plotter. Future versions of 
the postprocessing system will likely feature interactive color graphics which should 
greatly improve the usability of the modeling system as a research tool for studying 
atmospheric processes. Figs. 4-12 are examples of the output from three of the six 
postprocessing programs currently available. 
50 
M
R
SS
 2
.0
 
M
R
SS
 2
.0
 
50
0 
50
0 
T T 
VT
 
VT
 
21
00
 
21
00
 
G
M
T 
04
/0
2/
82
 
G
M
T 
04
/0
2/
82
 
IN
IT
 
IN
IT
 
12
00
 
12
00
 
G
M
T 
G
M
T 
04
/0
2/
82
 
04
/0
2/
82
 
FI
G
. 
4 
KA
SS
 FO
R
EC
AS
TS
 5
00
 M
E 
TE
M
PE
R
AT
U
R
ES
 (
D
EG
R
EE
S 
C
EL
C
IU
S)
 
5 
M
R
SS
 2
.0
 
50
0 
z 
VT
 
21
00
 
G
M
T 
04
/0
2/
82
 
IN
IT
 
12
00
 
G
M
T 
04
/0
2/
82
 
FI
G
, 
5 
M
AS
S F
O
R
EC
AS
TE
D
 5
00
 M
B 
H
EI
G
H
TS
 
(M
ET
ER
S)
 
AN
D
 
VO
R
TI
C
IT
Y 
(P
ER
 
SE
C
O
N
D
) 
0 
M
R
SS
 2
.0
 
50
0 
H
N
O
 
VT
 
21
00
 
G
M
T 
04
/0
2/
82
 
IN
IT
 
12
00
 
G
M
T 
04
/0
2/
82
 
FI
G
, 
6 
M
AS
S F
O
R
EC
AS
TE
D
 
50
0 
M
B 
W
IN
D
 
VE
C
TO
R
S 
AN
D
 
IS
O
TA
C
H
S 
(M
ET
ER
S 
PE
R
 
SE
C
O
N
D
) 
M
R
SS
 2
.0
 
50
0 
0M
E 
VT
 
21
00
 G
M
T 
04
/0
2/
82
 
IN
IT
 
12
00
 G
M
T 
04
/0
2/
82
 
FI
G
, 
7 
M
AS
S F
O
R
EC
AS
TE
D
 
50
0 
M
B 
VE
R
TI
C
AL
 
VE
LO
C
IT
Y 
(M
IC
R
O
BA
R
S 
PE
R
 
SE
C
O
N
D
) 
I 
5 
9 
13
 
17
 
21
 
25
 
29
 
33
 
37
 
41
 
45
 
99
 
53
 
57
 
61
 
65
 
69
 
73
 
77
 
61
 
85
 
39
 
93
 
97
lO
ll~
lO
91
l3
l1
71
2l
l~
 
FI
G
, 
8 
VE
R
TI
C
AL
 
C
R
O
SS
-S
EC
TI
O
N
 
LO
C
AT
O
R
 
M
AP
 
1 
I 
1)
 
17
 
21
 
25
 
26
 
33
 
37
 
41
 
I 
W
 
St
4 
17
 
61
 
#!
I 
89
 
73
 
~7
 
61
 
63
 
30
 
@
!j 
~7
10
11
tS
10
11
l3
11
71
21
12
5 
89
 
W
 
85
 
01
 
77
 
73
 
W
 
85
 
I1
 
57
 
53
 
¶9
 
45
 
41
 
87
 
99
 
28
 
25
 
21
 
17
 
13
 8 5 1
 
13
.5
 
13
.0
 
12
.5
 
12
.0
 
11
.5
 
11
.0
 
10
.5
 
10
.0
 
9.
5 
9.
0 
6.
5 
6.
0 
7.
5 
g 
‘,:
I 
6.
0 
5.
5 
5.
0 
4.
5 
4.
0 
3.
5 
3.
0 
2.
5 
2.
0 
1.
S 
1.
0 .S
 
.O
 
Q
H
 
- 
VT
 
21
00
 G
M
T 
04
/0
2/
82
 
IN
IT
 
12
00
 G
M
T 
04
/0
2/
82
 1
3.
5 
13
.0
 
12
.5
 
12
.0
 
11
.5
 
11
.0
 
10
.5
 
10
.0
 
9.
5 
9.
0 
6.
5 
8.
0 
7.
5 
7.
0 
6.
5 
2 
6.
0 
5.
5 
5.
0 
4.
5 
4.
0 
3.
5 
3.
0 
2.
5 
2.
0 
1.
5 
1.
0 .S
 
.O
 
FI
G
. 
9 
M
AS
S F
O
R
EC
AS
TE
D
 
VE
R
TI
C
AL
 
C
R
O
SS
-S
EC
TI
O
N
 
O
F 
R
EL
AT
IV
E 
H
U
M
ID
IT
Y 
(P
ER
C
EN
T)
 
13
.5
 
13
.0
 
12
.5
 
12
.0
 
1l
.S
 
11
.0
 
10
.5
 
10
.0
 
9.
5 
9.
0 
8.
5 
8.
0 
7.
5 
7.
0 
g 
6.
5 
6.
0 
5.
5 
5.
0 
4.
5 
4.
0 
3.
5 
3.
0 
2.
5 
2.
0 L.
5 
1.
0 .5
 
.O
 
13
.5
 
- 
13
.0
 
- 
12
.5
 
- 
12
.0
 
- 
11
.5
 
- 
11
.0
 
- 
10
.5
 
- 
10
.0
 
- 
9.
5 
- 
9.
0 
- 
8.
5 
- 
8.
0 
- 
7.
5 
- 
7.
0 
- 
6.
5 
2 
- 
6.
0 
- 
5.
5 
- 
5.
0 
- 
4.
5 
- 
4.
0 
- 
3.
5 
- 
3.
0 
- 
2,
s 
- 
2.
0 
- 
1.
5 
- 
1.
0 
J 
I- 
I 
I 
I 
I 
I 
i 
1 
I 
.5
 
.O
 
vv
 
VT
 
21
00
 
C
M
T 
04
/0
2/
82
 
IN
Il 
12
00
 
G
M
T 
04
/0
2/
U
 
FI
G
, 
10
 
Fl
AS
S F
O
R
EC
AS
TE
D
 
VE
R
TI
C
AL
 
C
R
O
SS
-S
EC
TI
O
N
 
O
F 
VE
R
TI
C
AL
 
VE
LO
C
IT
Y 
(M
 I 
C
R
O
BA
R
S 
PE
R
 
SE
C
O
N
D
) 
!5
0U
N
D
lN
t 
LE
C
AT
I0
N
S 
FI
G
, 
11
 
SO
U
N
D
IN
G
 
LO
C
AT
O
R
 M
AP
 
PT
 
N
U
M
BE
R
 
1.
 
ST
= 
PS
B 
cn
 
W
 
60
0 
VT
 
21
00
 
G
M
T 
04
/0
2/
82
 
IN
IT
 
12
00
 
G
M
T 
04
/0
2/
82
 
II~
~,
‘~
!‘~
I~
!‘~
i!‘
,!I
!“~
‘~
~’
 
I:.
‘: 
1.
‘. 
I’!
!““
:‘l
’!!
:!‘
!!!
l!.
“!!
‘“’
 
.,:
i!!
iq
 
16
15
1.
7 
XR
FR
C
E 
rn
Es
sL
M
 
96
8.
6 
L1
r-m
 
In
3 11
.7
 
m
m
! 
Tc
w
lrl
 
54
.2
 
LN
 
1 
TE
PP
 
IF
I 
47
.1
 
.7
0.
0 
-6
0.
0 
-5
0.
0 
-4
’2
.0
 
-3
3.
0 
-2
0.
0 
-!c
I.?
 
0.
0 
10
.0
 
20
.0
 
30
.0
 
40
.0
 
FI
G
, 
12
 
M
AS
S F
O
R
EC
AS
TE
D
 
SO
U
N
D
IN
G
 
12. REFERENCES 
Kaplan, M.L., J.W. Zack, V.C. Wong, and J.J. Tuccillo, 1982: Initial Results from a 
Mesoscale Atmospheric Simulation System and Comparisons with the AVE-SESAME I 
Data Set. Mon. Wea. Rev., 110, 1564-1590. 
Tuccillo, J.J., 1983: The Application of Pressure Gradient Force Averaging to the Euler- 
Backward Scheme. M.S. Thesis, Old Dominion University. 
60 


Memory efficient solution of the primitive equations for numerical weather prediction on the CYBER 205

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server