Circumventing Storage Limitations in Variational Data Assimilation Studies.
An application of Pontryagin's maximum principle data assimilation is used to blend possibly incomplete or nonuniformly distributed spatio-temporal observational data into geophysical models. Used extensively in engineering control theory applications, data assimilation has been introduced relatively recently into meteorological forecasting, natural-resource recovery modeling, and climate dynamics. Variational data assimilation is a promising assimilation technique in which it is assumed that the optimal state of the system is an extrema of a carefully chosen cost function. Provided that an adjoint model is available, the required model gradient can be computed by integrating the model forward and its adjoint backward. The gradient is then used to extremize the cost function with a suitable iterative or conjugate gradient solver.The problem addressed in this study is the explosive growth in both on-line computer memory and remote storage requirements of computing the gradient by the forward/adjoint technique which characterizes large-scale assimilation studies. Storage limitations impose severe limitation on the size of assimilation studies, even on the largest  computers. By using a recursive strategy, a schedule can be constructed that enables the forward/adjoint model runs to be performed in such a way that storage requirements can be traded for longer computational times.  This generally applicable strategy enables data assimilation studies on significantly larger domains than would otherwise be possible, given the particular hardware constraints, without compromising the outcome in any way. Furthermore, it is  shown that this tradeoff is indeed viable and that when the schedule is optimized, the storage and computational times grow at most logarithmically.
Introduction
Data assimilation has relatively recently become an important tool in many areas of
geophysics, such as weather and climate forecasting[1-6]; model sensitivity-to-parameter
studies[7,8], in the inclusion of field data sets into theoretical model-studies [9-11]. In
weather forecasting field data that may be spatially and/or temporally heteregeneous is
continuously blended into dynamical models as soon as the field data is available. the result
has been a significant improvement of the predictive capabilities of today's weather models
[1,12]. Ocean forecasting has thus far not met the comparable success of its atmospheric
counterpart. Reasons for this are (1) the spatial and temporal scales of the relevant oceanic
dynamics are several orders of magnitude smaller and larger, respectively; (2) oceanic data
gathering is at present very limited in coverage and sometimes of incompatible quality, (3)
boundary fluxes at the air/sea interface are poorly understood and yet have a major
influence on oceanic flows, and lastly, (4) the computing demands of oceanic forecasting
have only recently become marginally suitable for some but not all of the types of studies
at reasonable resolutions.
A specific approach to data assimilation is called variational data assimilation [12]. An
objective function is defined which provides a norm of the distance or misfit of the state
set to observational data. The state set may comprise model predictions, parameters,
boundary data and/or initial conditions. The misfit is usually weighted in order to account
for measurement errors, model uncertainties, etc. The object is to find the state set that
extremizes the objective function. This is usually done as a constrained optimization
problem which is generally solved iteratively by some extention of Newton's method or a
descent algorithm.
The optimization problem requires the computation of the gradient of the model with
respect to the state set. One of the alternative strategies that accomplishes the calculation
of the gradient is the "adjoint method" [3]. Provided an adjoint to the tangent linear
model exists, the process of computing the gradient involves integrating the original model
forward in time (the forward problem) recording the model's history, and then using the
history in the adjoint model to integrate backwards in time back to the point of origin
(the adjoint problem). Along the way the partial differentials that comprise the gradient
of the results at some t final with respect to the state set at some particular time step are
multiplied in reverse order until the adjoint model reaches the origin once again. By the
chain rule, the multiplication will yield the gradient, and it will do so at a computational
cost roughly twice that of the forward problem.
As described above, the adjoint method is what we will call the "conventional approach."
Its main advantage is its low computational cost. However, its disadvantage is that it
quickly encounters computer memory storage problems even in low resolution studies. In
this paper we present an alternative to the conventional approach which will circumvent in
a significant way the storage problems of the adjoint method at the expense of a possibly
greater, but manageable computational expense. The problem is motivated in Section
2. The alternative gradient method is presented in Section 3 and is compared to the
conventional approach. Section 4 demonstrates how such alternative is implemented in
practice in an ocean climate problem and we describe how it compares to the conventional
approach in terms of computational effort and memory usage. The last section summarizes
the findings, provides details of the strategy's computer implementation and details of
where to obtain code that implements the method.
[2]. Statement of the Problem
For the sake of clarity we will assume that the physical problem in question can be
modelled by an evolutionary equation. The physical m-dimensional real domain is R
with boundary @R. The discrete evolution parameter is so that physical
time marches according to i\Deltat. The goal is to find the state set which minimizes the
misfit of equation-to-field-data. The semi-discretized "forward problem" is
where the completely or partially unknown U and V are the initial and boundary data for
the state set which minimize an objective function. The "reverse problem" is the adjoint
of (2.1),
If the forward problem is a semi-discretization of an evolution equation, we think of u i and
with domain R \Theta T as vectors of the state variables and their adjoints, where i is the
time step index.
Equations (2.1) and (2.3) will be solved in some high level computer language such as
Fortran or C. Define
k as the set of computer memory addresses
required to represent the vector set fug and fu   g at index location i, so that u j and u
have temporary memory locations s j and s
k respectively. It is assumed that s
s
We will call this temporary computer storage medium as the
"register".
Let f and f   be the representations of F and F   respectively in some high level computer
program, or "program" for short. These take the form of subroutines, functions, etc. The
action of f . Define the m- and t- norms as the memory and
time of execution of some program Q as kQkm and kQk t respectively. As will be evident in
what follows, these norms will amount to simple direct sums. The register memory of the
state set is it will be safe to assume that kS   km - R. The other type of
memory that will play an important role in the analysis is the available memory external
to the program. This is usually some external storage device such as a memory disk, tape,
etc. For simplicity we shall call this recording device the "tape" and it will be assumed to
have a fixed memory size T . The distinction between a non-recording program procedure
f i and the same procedure which records the state set on tape will be indicated as -
f i . It
will be convenient to define the following specific m- and t- norms:
0-i-n
0-i-n
4 J.M. RESTREPO, G. K. LEAF, A. GRIEWANK
the maximum memory required to restore f i given S and the maximum computing time
(wall-clock time) to execute f i . It is worth noting that - is essentially fixed, regardless
of the number of processors while - can vary significantly depending on the number of
processors. Since f   is a linear mapping on S   , it can be assumed that
where -   and - refer respectively to analogous norms to (2.4) of f   and -
f , and the c's are
positive multiplicative constants. Note that kf i k - R since the subroutines may require
working registers.
In the discretization and coding of a typical evolution equation, for example of a climate
or metereology problem, we can identify f i as the collection of subroutines and functions
that take the state set from time t i to t i+1 (forward integration) in which kf i km and kf i k t
are approximately the same for each level thus equal to - and - respectively.
In the same fashion f
i is the collection of subroutines that take the state set from time t i to
(reverse integration) in which kf
are approximately the same for each
level thus equal to -   and -   respectively. Let us consider the memory and
the time norms of two strategies which may be used in the n\Gammastep gradient computation
by the adjoint method.
In one strategy the minimal memory norm is achieved by recording nothing on tape. It
requires stepping forward from u 0 to un using f i , followed by a single reverse step from un
to using f
n . Start again from u 0 forward to un\Gamma1 using f i followed by a reverse f
from un\Gamma1 to un\Gamma2 . The process is repeated until the reverse integration reaches step 0
once again. The t- and m- norms for this strategy are respectively
where only register memory is used. For simplicity we are ignoring here, as we will do
from now on, the register memory that is used for working arrays, etc. For an explicit
fourth-order Runge-Kutta time integration scheme for example, this register memory can
be significant, but can be easily accounted in the estimates provided.
Another strategy is the conventional approach, which steps forward from u 0 to un using
steps in reverse using f
reading the appropriate state variables from tape. The
time and memory norms for the latter strategy are
Hence the conventional approach yields the adjoint as a fixed multiple
of the time for the forward program. However, the tape grows linearly in both number
of steps and size of the state set, which for typical geophysical applications will quickly
overwhelm even the largest storage capabilities of computer facilities [13].
[3]. Recursive Adjoint Method
The recursive strategy is specifically designed to circumvent the storage limitations of
the conventional adjoint method at the expense of a larger computational effort. The
computational effort will be defined more precisely below, but for now it suffices to know
that the computational effort is directly proportional to the wall-clock time, which in turn
depends on the number of processors. One strategy that reduces the tape size is to produce
the gradient using the usual forward/adjoint sweep but recording less often than is really
required. While this alternative saves some tape space, it produces a degraded gradient. It
will be shown below that the gradient produced by the recursive method will be identical
to its nondegraded counterpart obtained in the conventional way.
The description that follows will present a heuristic explanation of the theoretical development
that appears in [14]. The basis of this strategy is to limit the tape size to
dR, , where d - n snapshots (snaps, for short) of states fug at any given point during
the program execution. This is done by carefully over-recording. It requires at most an
additional r-fold increase in additional full forward unrecorded computations or "reps".
The recursive strategy or "schedule" is not unique. However, from Theorem 6.1 due to
Griewank [14], among the partitioning algorithms the "binomial partitioning" schedule is
optimal. The theorem states that an n-step gradient calculation using the adjoint method
can be solved recursively using up to d - 0 snaps and at most r - 0 reps if and only if
d!r!
Note that n(d; d) and n(0; To illustrate the sense in which this
method is optimal we appeal to Stirling's formula, and find that for a fixed d or r,
or
To see more clearly the relationship between n and the number of snaps and reps, a
contour plot of ln n as a function of the number of snaps and reps based on (3.1) appears in

Figure

1. Since the values that the binomial takes are discrete, the contours appear jagged.
The figure clearly illustrates the logarithmic rate of growth given in (3.2). In particular,
when then they are growing as log 4 n.
The schedule for appears in Figure 2 and is worth
explaining in some detail. Note that
s
Along the horizontal is the number of reps and along the vertical the time step i. The
tree structure of the schedule is evident. Horizontal lines are drawn at locations in which
recordings are performed. As is evident, when reading the figure from left to right, there
6 J.M. RESTREPO, G. K. LEAF, A. GRIEWANK
reps3.054.577.6210.713.716.819.822.925.9
27.4 29 30.5

Figure

1. Contours of ln n versus snaps d and reps r.
are 5 self-similar groups or pennants. The top pennant and the first to be executed has
3 snaps at recording occurs at
s
and the recording at time step Execution requires a forward sweep from
56. The state at 51 is restored once more and a forward sweep to 55 follows. A
forward/adjoint from 55 to 56 and back again to 55 is executed then. The first pennant
is completely swept by repeating the last two steps until the adjoint reaches 51. State
36 is then restored and a forward sweep follows, recording at
the second uppper-most triangle is swept through, state 36 is recovered, a forward sweep
follows, recording at completing the first pennant state 1 is
restored, a forward sweep is initiated that ends at
s
recording along the way
at
s
+1. At this point, the schedule should be obvious. the last pennant is performed
when
s
2. Note that at no instant will the depth of the recording be more
than 3 records long, and in addition, that if the tape is thought of as a stack, the order of
the records is maintained due to its last-in-first-out nature. It is evident from the figure
that there are a total of 1 forward recorded sweep, 1 adjoint reverse sweep and r forward
unrecorded sweeps.
From

Figure

it may be concluded that the t-norm and m-norm of the recursive schedule
'i
'i
'i 3
'i
'i 5
'i 6
'i
'i 8
'i 9
'i
'i 11
'i
'i 12
'i
'i 14
'i 15
'i
'i 17
'i
'i 19
'i
'i 21
'i
'i 22
'i
'i
'i
'i
'i
'i 28
'i 29
'i
'i
'i
'i 33
'i 34
'i
'i 36
'i
'i 37
'i
'i
'i
'i 42
'i 43
'i 44
'i
'i
'i
'i 48
'i
'i
'i 52
'i 53
'i 54
'i
'i 56

Figure

2. Schedule for
are, respectively,
since dR. The first expression on the right hand side of equations (3.2) and (3.3) hold
generally for any n(d; r), d - 0 and r - 0 recursive adjoint problem and the far right hand
side for any general recursive adjoint problem involving the evolution equation typically
encountered in climate or metereology studies. Also note that if the number of reps r and
sweeps d are similar that
R
log 4 n:
8 J.M. RESTREPO, G. K. LEAF, A. GRIEWANK

Table

1. Schedule details for several sets of d, r and n.
Comparison of (3.2) and (2.6) lead to a working measure of the "computational effort"
which is directly proportional to the wall-clock time: a convenient measure is the total
number of forward steps. We shall employ this measure in this and in the following section,
in which a comparison between the recursive and the conventional approach is effected.
The performance of the recursive method compared to the conventional one may be
assessed graphically. Figure 3 illustrates the relation of the memory, measured in snaps, and
the wall-clock time, assuming it is proportional to the effort. The conventional approach
is represented by the left-most curve. All other curves represent different snap and rep
combinations. In both the conventional as well as the recursive case the memory required
to solve a problem will be equal to d- , where - is defined as before and depends on the
resolution and the number of spatial dimensions in the problem. The effort, on the other
hand, for the conventional case is basically n, while in the recursive strategy it depends on
the choice of snaps and reps. From left to right the recursive strategy curves correspond to
decreasing the number of snaps. The line-connected curve in the lower corner corresponds
to the case of snaps and reps being equal. The conventional case is in effect the limit of
snaps d equal to n in the recursive strategy. As can be surmized, the curves reflect the
previously-mentioned characteristic of the recursive method, that the effort increases for
the recursive method when fewer snaps are used. Hence, in practice, the user wants to
maximize the number of snaps in the calculation rather than the number of reps.
log(effort)41220283644snaps

Figure

3. Conventional versus Recursive strategy comparison. The
added effort due to increased reps r. From left to right, the conventional
case, then 3. The curve represented by stars
corresponds to r = d.

Figure

4 shows a comparison of the conventional strategy, the left-most solid curve, to
the recursive strategy with regards to the effort given by n. The slope of the recursive
curves gets closer to the slope of the conventional case the more snaps are used. Note
that in the conventional case the number of snaps is equal to n. Hence, this figure shows
the clear advantage of the recursive method with regards to memory. Specifically, whereas
an increase in n in the conventional case leads to an increase in tape usage, the recursive
strategy enables the user to consider a wider range of n for a fixed tape size dR. The
J.M. RESTREPO, G. K. LEAF, A. GRIEWANK
feasability of this latter strategy is dictated by the speed of the machine or the willingness
to pay for the higher effort involved. Compare this to the previous figure which shows the
price payed in higher wall-clock times, the smaller the number of snaps is employed. It may
be that the effort required in large problems is significant, but this must be weighed against
the fact that these problems may be simply impossible to consider with the conventional
strategy.
log(effort)2.06.010.0log(n)
conventional
d=2
d=3

Figure

4. Comparison of the conventional and recursive strategy.
Effort versus n. The memory requirement of the conventional case is
n. The recursive curves are labelled according to the number of snaps d
used. Natural logarithms are used.
[4]. Application to a Quasi-Geostrophic Ocean Problem
The recursive procedure's viability will be demonstrated by applying it to a quasi-geostrophic
model [15] [16], which was considered in Tziperman and Thacker's study [13],
hereon referred to as T&T. The dimensionless equations over a unit-square box in x and
y are
where /(x; are the streamfunction and the vorticity, - (x; y) is the wind
stress. J(\Delta; \Delta) is the symbol for the Jacobian of its arguments, and \Delta is the Laplacian
operator. The dimensionless real parameters R, ffl b , and ffl h are the Rossby number, the
bottom friction factor and the horizontal friction factor respectively. The state variables
evolve in time t and are subject to no-flux and no-stress boundary conditions at the edges
of the box.
The equations were discretized using multi-grid finite-difference techniques. The formulation
was first-order accurate in time and second-order accurate in space. In what follows
it will be understood that the state variables are only defined on the uniformly discretized
grid in x and y. We will omit explicit mention that these quantities are discretized in
space for the sake of clarity. On a discrete time grid i\Deltat, the state variables i i and / i
evolve to a steady state e
/. Following T&T, an assimilation problem is defined as
follows. The observational data will be the steady-state vorticity e
i which is independent
of time. The state set is taken to be the forcing term curl- , the initial vorticity i 0 , and
the parameters ffl b and ffl h . The observations e
are determined from a particular (fixed)
choice of friction factors effl b and effl h , initial vorticity i 0 and forcing curle- . The system is
then integrated forward in time until a steady-state is reached at which point the observations
are recorded. For purposes of this artificial assimilation problem, we now "forget"
the state set values which produced the observations. The task of the assimilation will
then be to reconstruct the state set which generated the observations. To this end, a cost
function is chosen which measure the fit of the model result to the observations. Since the
observations represent the steady-state, the cost function should measure the departure of
the model from steady-state as well as the departure from the observations. In T&T the
authors use the following discrete cost function:
where the sum is meant to indicate a sum over all the discrete values of the variables over
the unit box. The first term measures the deviation from the observations, while the second
term in conjunction with the first measures the deviation from steady-state. The matrices
C and D are the inverse of the covariance matrices of the observations. The final time
step n is arbitrary in this problem. It is chosen to be sufficiently large so that steady-state
is achieved. A small value of n reduces the computational cost per optimization iteration
however, it increases the number of optimization iterations. Since the number of recorded
histories depends on the number of time steps n, the storage requirements are reduced
when n is small. In fact, in T&T In their experiment such a choice is possible since
the assimilation occurs at just one time level. In the general case, assimilations may occur
at multiple time levels, in which case the number of time steps used is determined by the
problem and cannot be arbitrarily chosen.
The optimization task is to find the state set fcurl-; for which H n is a minimum
subject to the constraints of the model equations. A common strategy for computing
the minimum is to introduce Lagrange multipliers and the corresponding Lagrange functions
for which we seek an unconstrained extremum. A gradient-based iterative algorithm
such as the conjugate gradient method was then applied to this unconstrained problem.
J.M. RESTREPO, G. K. LEAF, A. GRIEWANK
For the discrete quasi-geostrophic model, the Lagrange function has the form
ae @i
@t
@x
oe
The descent algorithm requires the calculation of the gradient of L n with respect to the
state set. The gradient involve the Lagrange multipliers f- which are determined from
the gradients of L n with respect to fi g. Equating these gradients to zero generates
the adjoint equations for f- which may be symbolically expressed as
where \Psi is the forcing term arising from the gradients of the cost function with respect
to g. The discrete adjoint equations are integrated backwards in time to generate
the Lagrange multipliers - i used in computing the gradients of the cost function as needed
by the conjugate gradient procedure. Thus, in the conventional approach, each conjugate
gradient iteration requires a forward integration of n steps, which generates the value
of the cost function, followed by a backward integration of the adjoint equations. This
adjoint integration generates the gradients used in the conjugate gradient iteration. It is
noted that the state set is required to effect the calculation of the Lagrange multipliers.
Ordinarily the state set at each time step i is required for the calculation in the adjoint
problem, hence these are recorded in the forward sweep. Since only the state variables are
time-dependent in this particular problem, we need only record the state variables i
at each time step. The remaining components of the state set need to be recorded only
once during the forward-backwards sweep. The observations were synthesized by running
the discretized version of (4.1) to steady-state using
To demonstrate the performance of the recursive forward-backward integration strategy
for the calculation of the gradient we compared model runs of this experiment using T&T's
original multigrid Fortran code against a version of the code which was identical in all
respects to T&T's code except for a subroutine that generates the schedule and minor
modifications to the program to enable us to carry the schedule out. As a first step, we
thoroughly verified that our program results yielded identical results to the conventional
case. The wall-clock time was negligibly higher for the recursive program running in
the conventional mode, reflecting the additional computational expense of generating the
schedule.
In the experiments to be reported, the optimality tolerance for the NAG conjugate
gradient routine was set to 10 \Gamma3 in all model runs. The square of / summed over
the box was the residual in the calculation. The forward run used to create the observations
stepped in time until the residual was below 10 \Gamma7 . The multigrid depth was fixed at 4
levels for all experiments that follow. The codes were executed on a Sparc 10/51 running
SunOS 4.1.3U1. The Fortran Sun compiler used was Fortran Version 1.4 with optimization
memory (floats)100300500700n

Figure

5. Comparison of the conventional and recursive strategy on
the T&T problem. Tape versus n. In the recursive strategy the snaps
was held fixed at 5. The recursive strategy has a fixed tape length
of 1089 floats.
flags turned off. All runs were performed in double-precision arithmetic. Wall-clock times
reported encompass the solution to the full problem. In all experiments performed the
answers from both strategies were identical.

Figure

5 shows a comparison of tape usage for the conventional and the recursive strat-
egy. In the recursive trials the snap count was held fixed at 5, explaining why its curve
for tape usage is a vertical straight line. As mentioned previously, for the conventional
case the tape usage is proportional to the number of time steps n. From Figure 5 the tape
for the conventional case. It follows from this experiment that if there was a
fixed amount of tape on a particular machine, that the conventional approach would be
impossible for sufficiently large problems. Figure 6 shows the wall-clock time for the same
experiment. In all trials the conjugate gradient procedure converged in three iterations.
The conventional wall-clock time seconds, the recursive version took
14 J.M. RESTREPO, G. K. LEAF, A. GRIEWANK
time (seconds)100300500700n

Figure

6. Comparison of the conventional and recursive strategy for
the T&T problem. Wall-clock versus n. In the recursive strategy the
snaps was held fixed at 5. The rightmost curve corresponds to the
recursive strategy.
longer to complete and it's growth is not linear. Table 2 contains further information on
this particular set of trials.
Another typical problem in which the recursive strategy would prove useful is illustrated
by the following example. In this assimilation problem we have some flexibility in choosing
the integration time, since the only requirement is that it must be longer than n   , where
is the minimum number of steps for a steady-state solution. This is precisely the type of
problem investigated in Tziperman and Thacker's study [13]. Suppose that for a particular
resolution the problem "fits" and thus can be solved on a particular machine using the
conventional approach. In order to double the spatial resolution, the conventional strategy
would require a four-fold increase in tape storage. The doubly-resolved experiment no
longer could be performed on this particular machine. However, the problem could be
solved using the recursive approach so long as the maximum tape length is not exceeded.

Table

2. Ratio of the wall-clock time for the
recursive conventional approach versus
n and number of reps for the T&T problem.
Suppose that the maximum tape length on this machines is 60984 floats. This is the actual
requirement of the singly-resolved T&T problem with spatial grid
with 4 refinement levels. Table 3 provides the results of several runs using the recursive
strategy for the doubly-resolved problem. Supposing that the conventional procedure could
be could be carried out, for the tape length for the doubly-resolved problem would
be 236600 floats and it would have taken 40:42 seconds to execute. The table demonstrates
that the doubly-resolved problem can be succesfully carried out in two to four times the
amount of time that it would take to run the conventional procedure assuming that it
could be possible to compute conventionally in the first place.

Table

3. Wall-clock time and tape length for
the recursive conventional approach versus n in
the T&T problem for a doubling of resolution.
time (secs) tape (floats) snaps reps
[4]. Conclusions
We have shown in this study how a recursive strategy due to Griewank [14] for the
adjoint-method calculation of the gradient may be applied to large-scale geophysical variational
data assimilation studies. The main result is that significantly larger assimilation
studies can be performed using this recursive strategy than is possible using the conventional
forward-adjoint methods, given the physical limitations of available computer storage
hardware, at the expense of a reasonable additional computational effort (or wall-clock
time). Furthermore, the procedure yields the gradient with no degradation as compared
J.M. RESTREPO, G. K. LEAF, A. GRIEWANK
to the conventional approach.
In theory, when the number of snaps and reps, i.e. the number of storage units (measured
in R), and the number of additional unrecorded forward runs, are equal, these are both
bounded by log 4 n, where n is the number of time steps in the evolution equation. In
practice, the strategy is best used by picking the maximum number of snaps that the
particular computer hardware can manage, thus minimizing the number of reps.
It is remarked that
s
are elements of Z + , is not a
complete covering of Z . Furthermore, the distance between neighboring elements of the
set gets larger the larger n is. The reader may thus be wondering what is the strategy
when
options may be viable: (1) pick an n   2 Z n which is greater than
or equal to the required number of time steps in the problem, fixing d, say, and adjusting
r so that n   is slightly greater than n; (2) pick a finer (or coarser) time discretization;
Circumventing storage limitations by far outweigh the relatively
minor complication presented by the fact that
Option (1) is the simplest and is
suitable in the T&T problem. We foresee that the recursive strategy would be simplest
to implement in problems like steady-state climate, hydrology, oil-recovery, etc., and least
straightforward in problems like metereological forecasting, earthquake forecasting, and
the like. For the forecasting problems perhaps an alternative partitioning algorithm with
a denser set Z n may be considered more appropriate.
Insofar as computer program design, the best strategy for large-scale problems is to construct
programs that are as compute-intensive as possible and the least memory-intensive.
This yields the greatest variation in the computational effort for any given choice of snaps
and reps. This is especially true in parallelized programs because the computational effort
will drop the greater the number of processors used, whereas the storage requirements
remain fixed independent of the number of processors.
The implementation of the recursive strategy requires minimal modifications of conventional
codes that compute forward and adjoint problems. The requirements are that four
modules be provided, namely (1) a forward module that runs without recording the state
set between a specified starting and an ending time step; (2) a module that computes a
single unrecorded forward and a single adjoint step, given a specific time step; (3) a module
that records to tape the state set at the current time step; and (4) a module that retrieves
from tape the last recorded state set. An additional module, which is to be considered the
driver, runs the above-mentioned modules according to the recursive schedule. The driver
requires as input the total number of time steps, the number of snaps and the number of
reps.
One approach in the implementation of the schedule driver is to have the schedule
computed only once at the top of the program. The schedule instructions are saved in
integer arrays which are then called in sequence to drive the four modules. The benefit of
precomputing the schedule is not warranted in some applications, since the schedule module
increases insignificantly the overall computational effort. The preferred alternative is to
use the schedule driver to control the above-mentioned modules thus not wasting register
memory for the schedule arrays needed in the first approach that could otherwise be used
in the adjoint problem. An estimate of the additional memory for the integer schedule
arrays of the first approach is as follows: a "schedule array" with the instruction directives
of size 2rn is required, plus one or two arrays of similar size which direct the recording
and accessing of snaps from tape. The total register overhead is then in the order of 4rn
integers. The user's particular application will clearly dictate which of these alternatives
works best.
This schedule driver is available via anonymous ftp from info.mcs.anl.gov. The file is
called /pub/tech reports/restrepo/schedule.tar.Z. Alternatively, the schedule software
is available in either Fortran or C versions from the word-wide-web in the software
section of http://www.mcs.anl.gov/people/restrepo/index.html.
J.M. RESTREPO, G. K. LEAF, A. GRIEWANK



--R


Metereological Data Assimilation for Oceanographers.


Variational Four Dimensional Analysis Using Quasi-Geostrophic Con- straints
Finding the Steady State of a General Circulation Model Through Data Assimilation: Application to the North Atlantic Ocean
Newton's Method with a Model Trust Region Modifications
Solution of Nonlinear Finite Difference Ocean Models by Optimization Methods with Sensitivity and Observational Strategy Analysis
The Role of Integration Time in Determining a Steady State through Data Assimilation
An Adjoint Method for Obtaining the Most Rapidly Growing Perturbation to Oceanic Flows
Fitting Dynamics to Data
Variational Assimilation of Meteorological Observations with the Adjoint Vorticity Equation
An Optimal Control-Adjoint Equations Approach to Studying the Oceanic General Circulation
Achieving Logarithmic Growth in Temporal and Spatial Complexity in Reverse Automatic Differentiation
Driven Ocean Circulation- Part 2

--TR

--CTR
Andreas Griewank , Andrea Walther, Algorithm 799: revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation, ACM Transactions on Mathematical Software (TOMS), v.26 n.1, p.19-45, March 2000
Patrick Heimbach , Chris Hill , Ralf Giering, An efficient exact adjoint of the parallel MIT general circulation model, generated via automatic differentiation, Future Generation Computer Systems, v.21 n.8, p.1356-1371, October 2005
Yanhua Cao , Jiang Zhu , Zhendong Luo , I. M. Navon, Reduced-Order Modeling of the Upper Tropical Pacific Ocean Model using Proper Orthogonal Decomposition, Computers & Mathematics with Applications, v.52 n.8-9, p.1373-1386, October, 2006
Matthias Heinkenschloss, A time-domain decomposition iterative method for the solution of distributed linear quadratic optimal control problems, Journal of Computational and Applied Mathematics, v.173 n.1, p.169-198, 1 January 2005
gradient;variational data assimilation;storage;recursion;adjoint model;natural resources recovery;climate;meteorology