1 |
% $Header: /u/gcmpack/manual/s_examples/global_oce_optim/global_oce_estimation.tex,v 1.16 2010/08/27 13:25:32 jmc Exp $ |
2 |
% $Name: $ |
3 |
|
4 |
\section[Global Ocean State Estimation Example]{Global Ocean State Estimation at $4^\circ$ Resolution} |
5 |
%\label{www:tutorials} |
6 |
\label{sec:eg-global_state_estimate} |
7 |
\begin{rawhtml} |
8 |
<!-- CMIREDIR:eg-global_state_estimate: --> |
9 |
\end{rawhtml} |
10 |
\begin{center} |
11 |
(in directory: {\it verification/tutorial\_global\_oce\_optim/}) |
12 |
\end{center} |
13 |
|
14 |
\subsection{Overview} |
15 |
|
16 |
This experiment illustrates the optimization capacity of the MITgcm: here, |
17 |
a high level description. |
18 |
|
19 |
In this tutorial, a very simple case is used to illustrate the optimization |
20 |
capacity of the MITgcm. Using an ocean configuration with realistic geography |
21 |
and bathymetry on a $4\times4^\circ$ spherical polar grid, we estimate a |
22 |
time-independent surface heat flux adjustment $Q_\mathrm{netm}$ that attempts |
23 |
to bring the model climatology into consistency with observations (Levitus |
24 |
dataset, \cite{lev:94a}). The files for this experiment can be found in the |
25 |
verification directory under tutorial\_global\_oce\_optim. |
26 |
|
27 |
This adjustment $Q_\mathrm{netm}$ (a 2D field only function of longitude and |
28 |
latitude) is the control variable of an optimization problem. It is inferred |
29 |
by an iterative procedure using an `adjoint technique' and a least-squares |
30 |
method (see, for example, \cite{stam-etal:02} and \cite{fer-eta:05}). |
31 |
|
32 |
The ocean model is run forward in time and the quality of the solution is |
33 |
determined by a cost function, $J_1$, a measure of the departure of the model |
34 |
climatology from observations: |
35 |
\begin{equation}\label{cost_temp} |
36 |
J_1=\frac{1}{N}\sum_{i=1}^N \left[ \frac{\overline{T}_i-\overline{T}_i^{lev}}{\sigma_i^T}\right]^2 |
37 |
\end{equation} |
38 |
where $\overline{T}_i$ and $\overline{T}_i^{lev}$ are, respectively, the model |
39 |
and observed potential temperature at each |
40 |
grid point $i$. The differences are weighted by an {\it a priori} uncertainty |
41 |
$\sigma_i^T$ on observations (as provided by \cite{lev:94a}). The error |
42 |
$\sigma_i^T$ is only a function of depth and varies from 0.5 at the surface to |
43 |
0.05~K at the bottom of the ocean, mainly reflecting the decreasing |
44 |
temperature variance with depth (Fig. \ref{Error}a). A value of $J_1$ of |
45 |
order 1 means that the model is, on average, within observational |
46 |
uncertainties. |
47 |
|
48 |
The cost function also places constraints on the adjustment to insure it is |
49 |
"reasonable", i.e. of order of the uncertainties on the observed surface heat |
50 |
flux: |
51 |
\begin{equation} |
52 |
J_2 = \frac{1}{N} \sum_{i=1}^N \left[\frac{Q_\mathrm{netm}}{\sigma^Q_i} \right]^2 |
53 |
\end{equation} |
54 |
where $\sigma^Q_i$ are the {\it a priori} errors on the observed heat flux as |
55 |
estimated by Stammer et al. (2002) from 30\% of local root-mean-square |
56 |
variability of the NCEP forcing field (Fig \ref{Error}b). |
57 |
|
58 |
The total cost function is defined as $J=\lambda_1 J_1+ \lambda_2 J_2$ where |
59 |
$\lambda_1$ and $\lambda_2$ are weights controlling the relative contribution |
60 |
of the two components. The adjoint model then yields the sensitivities |
61 |
$\partial J/\partial Q_\mathrm{netm}$ of $J$ relative to the 2D fields |
62 |
$Q_\mathrm{netm}$. Using a line-searching algorithm (\cite{gil-lem:89}), |
63 |
$Q_\mathrm{netm}$ is adjusted then in the sense to |
64 |
reduce $J$ --- the procedure is repeated until convergence. |
65 |
|
66 |
%The configuration is identical |
67 |
%to the ``Global Ocean circulation'' tutorial where more details can be found. |
68 |
|
69 |
Fig. \ref{Results} shows the results of such an optimization. The |
70 |
model is started from rest and from January-mean temperature and salinity |
71 |
initial conditions taken from the Levitus dataset. The experiment is run a year |
72 |
and the averaged temperature over the whole run (i.e. annual mean) is used |
73 |
in the cost function (\ref{cost_temp}) to evaluate the model\footnote{Because of |
74 |
the daily automatic testing, the experiment as found in the repository is set-up |
75 |
with a very small number of time-steps. To reproduce the results shown here, one |
76 |
needs to set nTimeSteps = 360 and lastinterval=31104000 (both corresponding to a |
77 |
year, see section \ref{sectRun} for further details).}. Only the |
78 |
top 2 levels are used. The first guess $Q_\mathrm{netm}$ is chosen to be |
79 |
zero. The weights $\lambda_1$ and $\lambda_2$ are set to 1 and 2, respectively. |
80 |
The total cost function converges after 15 iterations, decreasing from 6.1 to |
81 |
2.7 (the temperature contribution decreases from 6.1 to 1.8 while the heat |
82 |
flux one increases from 0 to 0.42). The right panels of Fig. (\ref{Results}) |
83 |
illustrate the evolution of the temperature error at the surface from |
84 |
iteration 0 to iteration 15. Unsurprisingly, the largest errors at iteration 0 |
85 |
(up to 6$^\circ$C, top left panels) are found in the Western boundary |
86 |
currents. After optimization, the departure of the model temperature from |
87 |
observations is reduced to 1$^\circ$C or less almost everywhere except in the |
88 |
Pacific Equatorial Cold Tongue. Comparison of the initial temperature |
89 |
error (top, right) and heat flux adjustment (bottom, left) shows that the |
90 |
system basically increased the heat flux out of the ocean where temperatures |
91 |
were too warm and vice-versa. Obviously, heat flux uncertainties are not the |
92 |
sole responsible for temperature errors and the heat flux adjustment partly |
93 |
compensates the poor representation of narrow currents (Western boundary |
94 |
currents, Equatorial currents) at $4\times4^\circ$ resolution. This is |
95 |
allowed by the large {\it a priori} error on the heat flux (Fig. \ref{Error}). |
96 |
The Pacific Cold Tongue is a counter example: there, heat fluxes uncertainties |
97 |
are fairly small (about 20~W.m$^2$), and a large temperature errors |
98 |
remains after optimization. |
99 |
|
100 |
In the following, section 2 describes in details the implementation of the |
101 |
control variable $Q_\mathrm{netm}$, the cost function $J$ and the I/O required |
102 |
for the communication between the model and the line-search. Instructions to |
103 |
compile the MITgcm and its adjoint and the line-search algorithm are given in |
104 |
section 3. The method used to run the experiment is described in section 4. |
105 |
|
106 |
\begin{figure} [tpb] |
107 |
\begin{center} |
108 |
\includegraphics[width=\textwidth,height=.3\textheight]{s_examples/global_oce_optim/Error.eps} |
109 |
\caption{{\it A priori} errors on potential temperature (left, in $^\circ$C) and |
110 |
surface heat flux (right, in W~m$^{-2}$) used to compute the cost |
111 |
terms $J_1$ and $J_2$, respectively.} |
112 |
\label{Error} |
113 |
\end{center} |
114 |
\end{figure} |
115 |
|
116 |
\begin{figure} [tpb] |
117 |
\begin{center} |
118 |
\includegraphics[width=\textwidth,height=.3\textheight]{s_examples/global_oce_optim/Tutorial_fig.eps} |
119 |
\caption{Initial annual mean surface heat flux (top right in W.m$^{-2}$) and |
120 |
adjustment obtained at iteration 15 (bottom right). Averaged difference |
121 |
between model and observed potential temperatures at the surface (in $^\circ$C) |
122 |
before optimization (iteration 0, top right) and after optimization |
123 |
(iteration 15, bottom right). Contour intervals for heat flux and temperature |
124 |
are 25~W.m$^{-2}$ and 1$^\circ$C, respectively. A positive flux is out of the |
125 |
ocean.} |
126 |
\label{Results} |
127 |
\end{center} |
128 |
\end{figure} |
129 |
|
130 |
\subsection{Implementation of the control variable and the cost function} |
131 |
|
132 |
One of the goal of this tutorial is to illustrate how to implement a new |
133 |
control variable. Most of this is fairly generic and is done in the ctrl |
134 |
and cost packages found in the pkg/ directory. The modifications can be |
135 |
tracked by the CPP option ALLOW\_HFLUXM\_CONTROL or the comment |
136 |
cHFLUXM\_CONTROL. The more specific modifications required for the experiment |
137 |
are found in verification/tutorial\_global\_oce\_optim/code\_ad. Here follows |
138 |
a brief description of the implementation. |
139 |
|
140 |
\subsubsection{The control variable} |
141 |
|
142 |
The adjustment $Q_\mathrm{netm}$ is activated by setting |
143 |
ALLOW\_HFLUXM\_CONTROL to "define" in ECCO\_OPTIONS.h. |
144 |
|
145 |
It is first implemented as a ``normal'' forcing variable. It is defined in |
146 |
FFIELDS.h, initialized to zero in ini\_forcing.F, and then used in |
147 |
external\_forcing\_surf.F. $Q_\mathrm{netm}$ is made a control variable in |
148 |
the ctrl package by modifying the following subroutines: |
149 |
|
150 |
\begin{itemize} |
151 |
\item ctrl\_init.F where $Q_\mathrm{netm}$ is defined as the control variable |
152 |
number 24, |
153 |
|
154 |
\item ctrl\_pack.F which writes, at the end of each iteration, the sensitivity |
155 |
of the cost function $\partial J/\partial Q_\mathrm{netm}$ in to a file to be |
156 |
used by the line-search algorithm, |
157 |
|
158 |
\item ctrl\_unpack.F which reads, at the start of each iteration, the updated |
159 |
adjustment as provided by the line-search algorithm, |
160 |
|
161 |
\item ctrl\_map\_forcing.F in which the updated adjustment is added to the |
162 |
first guess $Q_\mathrm{netm}$. |
163 |
\end{itemize} |
164 |
|
165 |
Note also some minor changes in ctrl.h, ctrl\_readparams.F, and ctrl\_dummy.h |
166 |
(xx\_hfluxm\_file, fname\_hfluxm, xx\_hfluxm\_dummy). |
167 |
|
168 |
\subsubsection{Cost functions} |
169 |
|
170 |
The cost functions are implemented using the {\it cost} package. |
171 |
|
172 |
\begin{itemize} |
173 |
|
174 |
\item The temperature cost function $J_1$ which measures the drift of the mean |
175 |
model temperature from the Levitus climatology is implemented in cost\_temp.F. |
176 |
It is activated by ALLOW\_COST\_TEMP in ECCO\_OPTIONS.h. It requires the mean |
177 |
temperature of the model which is obtained by accumulating the temperature in |
178 |
cost\_tile.F (called at each time step). |
179 |
The value of the cost function is stored in {\it objf\_temp} and its weight |
180 |
$\lambda_1$ in {\it mult\_temp}. |
181 |
|
182 |
\item The heat flux cost function, penalizing the departure of the surface |
183 |
heat flux from observations is implemented in cost\_hflux.F, and activated by |
184 |
the key ALLOW\_COST\_HFLUXM in ECCO\_OPTIONS.h. The value of the cost |
185 |
function is stored in {\it objf\_hfluxm} and its weight $\lambda_2$ in |
186 |
{\it mult\_hfluxm}. |
187 |
|
188 |
\item The subroutine cost\_final.F calls the cost\_functions subroutines |
189 |
and make the (weighted) sum of the various contributions. |
190 |
|
191 |
\item The various weights used in the cost functions are read in |
192 |
cost\_weights.F. The weight of the cost functions are read in |
193 |
cost\_readparams.F from the input file data.cost. |
194 |
|
195 |
\end{itemize} |
196 |
|
197 |
|
198 |
\subsection{Code Configuration} |
199 |
%\label{www:tutorials} |
200 |
\label{sec:eg_globest_code_config} |
201 |
|
202 |
The model configuration for this experiment resides under the directory |
203 |
{\it verification/tutorial\_global\_oce\_optim/}. The experiment files in |
204 |
code\_ad/ and input\_ad/ contain the code customizations and parameter |
205 |
settings. Most of them are identical to those used in the Global Ocean |
206 |
( experiment {\it tutorial\_global\_oce\_latlon}). Below, we describe some of |
207 |
the customizations required for this experiment. |
208 |
|
209 |
\subsubsection{Compilation-time customizations in {\it code\_ad/}} |
210 |
|
211 |
In ECCO\_CPPOPTIONS.h: |
212 |
|
213 |
\begin{itemize} |
214 |
\item define ALLOW\_ECCO\_OPTIMIZATION |
215 |
|
216 |
\item define ALLOW\_COST, ALLOW\_COST\_TEMP, and ALLOW\_COST\_HFLUXM |
217 |
|
218 |
\item define ALLOW\_HFLUXM\_CONTROL |
219 |
\end{itemize} |
220 |
|
221 |
\subsubsection{Running-time customizations in {\it input\_ad/}}\label{sectRun} |
222 |
|
223 |
\begin{itemize} |
224 |
|
225 |
\item {\it data}: note the smaller {\it cg2dTargetResidual} than in the |
226 |
forward-only experiment, |
227 |
|
228 |
\item {\it data.optim} specifies the iteration number, |
229 |
|
230 |
\item {\it data.ctrl} is used, in particular, to specify the |
231 |
name of the sensitivity and adjustment files associated to a control |
232 |
variable, |
233 |
|
234 |
\item {\it data.cost}: parameters of the cost functions, in particular |
235 |
{\it lastinterval} specifies the length of time-averaging for the model |
236 |
temperature to be used in the cost function (\ref{cost_temp}), |
237 |
|
238 |
\item {\it data.pkg}: note that the Gradient Check package is turned on by |
239 |
default (useGrdchk=.TRUE.), |
240 |
|
241 |
\item {\it Err\_hflux.bin} and {\it Err\_levitus\_15layer.bin} are the |
242 |
files containing the heat flux and potential temperature uncertainties, |
243 |
respectively. |
244 |
|
245 |
\end{itemize} |
246 |
|
247 |
\subsection{Compiling} |
248 |
|
249 |
The optimization experiment requires two executables: 1) the |
250 |
MITgcm and its adjoint ({\it mitgcmuv\_ad}) and 2) the line-search |
251 |
algorithm ({\it optim.x}). |
252 |
|
253 |
\subsubsection{Compilation of MITgcm and its adjoint: {\it mitcgmuv\_ad}} |
254 |
|
255 |
Before compiling, first note that, in the directory code\_ad/, two files |
256 |
must be updated: |
257 |
\begin{itemize} |
258 |
\item code\_ad\_diff.list which lists new subroutines to be compiled |
259 |
by the TAF software (cost\_temp.F and cost\_hflux.F here), |
260 |
|
261 |
\item the adjoint\_hfluxm files which provides a list of the control variables |
262 |
and the name of cost function to the TAF software. |
263 |
|
264 |
\end{itemize} |
265 |
|
266 |
Then, in the directory build\_ad/, type: |
267 |
\begin{verbatim} |
268 |
% ../../../tools/genmake2 -mods=../code\_ad -adof=../code\_ad/adjoint\_hfluxm |
269 |
% make depend |
270 |
% make adall |
271 |
\end{verbatim} |
272 |
to generate the MITgcm executable {\it mitgcmuv\_ad}. |
273 |
|
274 |
\subsubsection{Compilation of the line-search algorithm: {\it optim.x}} |
275 |
|
276 |
This is done from the directories lsopt/ and optim/ (under MITgcm/). In lsopt/, |
277 |
unzip the blash1 library adapted to your platform, and change the Makefile |
278 |
accordingly. Compile with: |
279 |
\begin{verbatim} |
280 |
% make all |
281 |
\end{verbatim} |
282 |
(more details in lsopt\_doc.txt) |
283 |
|
284 |
In optim/, the path of the directory where {\it mitgcm\_ad} was compiled |
285 |
must be specified in the Makefile in the variable INCLUDEDIRS. The file name |
286 |
of the control variable (xx\_hfluxm\_file here) must be added to the name list |
287 |
read by optim\_num.F. Then use |
288 |
\begin{verbatim} |
289 |
% make depend |
290 |
\end{verbatim} |
291 |
and |
292 |
\begin{verbatim} |
293 |
% make |
294 |
\end{verbatim} |
295 |
to generate the line-search executable {\it optim.x}. |
296 |
|
297 |
\subsection{Running the estimation} |
298 |
|
299 |
Copy the {\it mitgcmuv\_ad} executable to input\_ad/ and {\it optim.x} to the |
300 |
subdirectory input\_ad/OPTIM/. Move into input\_ad/. The first iteration is |
301 |
somewhat particular and is best done "by hand" while the following iterations |
302 |
can be run automatically (see below). Check that the iteration number is set |
303 |
to 0 in data.optim and run the MITgcm: |
304 |
\begin{verbatim} |
305 |
% ./mitgcmuv_ad |
306 |
\end{verbatim} |
307 |
|
308 |
The output files adxx\_hfluxm.0000000000.* and xx\_hfluxm.0000000000.* contain |
309 |
the sensitivity of the cost function to $Q_\mathrm{netm}$ and the adjustment |
310 |
to $Q_\mathrm{netm}$ (zero at the first iteration), respectively. Two other |
311 |
files called costhflux\_tut\_MITgcm.opt0000 and ctrlhflux\_tut\_MITgcm.opt0000 |
312 |
are also generated. They essentially contain the same information as the |
313 |
adxx\_.hfluxm* and xx\_hfluxm* files, but in a compressed format. These two files |
314 |
are the only ones involved in the communication between the adjoint model |
315 |
{\it mitgcmuv\_ad} and the line-search algorithm {\it optim.x}. Only at the first |
316 |
iteration, are they both generated by {\it mitgcmuv\_ad}. Subsenquently, |
317 |
costhflux\_tut\_MITgcm.opt$n$ is an output of the adjoint model at |
318 |
iteration $n$ and an input of the line-search. The latter returns an updated |
319 |
adjustment in ctrlhflux\_tut\_MITgcm.opt$n+1$ to be used as an input of the |
320 |
adjoint model at iteration n+1. |
321 |
|
322 |
At the first iteration, move costhflux\_tut\_MITgcm.opt0000 and |
323 |
ctrlhflux\_tut\_MITgcm.opt0000 to OPTIM/, move into this directory and link data.optim |
324 |
and data.ctrl locally: |
325 |
\begin{verbatim} |
326 |
% cd OPTIM/ |
327 |
% ln -s ../data.optim . |
328 |
% ln -s ../data.ctrl . |
329 |
\end{verbatim} |
330 |
The target cost function {\it fmin} needs to be specified too: as a rule of thumb, |
331 |
it should be about 0.95-0.90 times the value of the cost function at |
332 |
the first iteration. This value is only used at the first iteration and does |
333 |
not need to be updated afterward. However, it implicitly specifies the |
334 |
``pace'' at which the cost function is going down (if you are lucky and it does |
335 |
indeed diminish!). More in the ECCO section maybe ? |
336 |
|
337 |
Once this is done, run the line-search algorithm: |
338 |
\begin{verbatim} |
339 |
% ./optim.x |
340 |
\end{verbatim} |
341 |
which computes the updated adjustment for iteration 1, ctrlhflux\_tut\_MITgcm.opt0001. |
342 |
|
343 |
The following iterations can be executed automatically using the shell |
344 |
script {\it cycsh} found in input\_ad/. This script will take care of changing |
345 |
the iteration numbers in the data.optim, launch the adjoint model, clean and |
346 |
store the outputs, move the costhflux* and ctrlhflux* files, and run the |
347 |
line-search algorithm. Edit {\it cycsh} to specify the prefix of the |
348 |
directories used to store the outputs and the maximum number of iteration. |
349 |
|