s_examples/global_oce_optim/global_oce_estimation.tex

% $Header: /u/gcmpack/manual/s_examples/global_oce_optim/global_oce_estimation.tex,v 1.16 2010/08/27 13:25:32 jmc Exp $
% $Name:  $

\section[Global Ocean State Estimation Example]{Global Ocean State Estimation at $4^\circ$ Resolution}
%\label{www:tutorials}
\label{sec:eg-global_state_estimate}
\begin{rawhtml}
<!-- CMIREDIR:eg-global_state_estimate: -->
\end{rawhtml}
\begin{center}
(in directory: {\it verification/tutorial\_global\_oce\_optim/})
\end{center}

\subsection{Overview}

This experiment illustrates the optimization capacity of the MITgcm: here,
a high level description.

In this tutorial, a very simple case is used to illustrate the optimization
capacity of the MITgcm. Using an ocean configuration with realistic geography
and bathymetry on a $4\times4^\circ$ spherical polar grid, we estimate a
time-independent surface heat flux adjustment $Q_\mathrm{netm}$ that attempts
to bring the model climatology into consistency with observations (Levitus
dataset, \cite{lev:94a}). The files for this experiment can be found in the
verification directory under tutorial\_global\_oce\_optim.

This adjustment $Q_\mathrm{netm}$ (a 2D field only function of longitude and
latitude) is the control variable of an optimization problem. It is inferred
by an iterative procedure using an `adjoint technique' and a least-squares
method (see, for example, \cite{stam-etal:02} and \cite{fer-eta:05}).

The ocean model is run forward in time and the quality of the solution is
determined by a cost function, $J_1$, a measure of the departure of the model
climatology from observations:
\begin{equation}\label{cost_temp}
J_1=\frac{1}{N}\sum_{i=1}^N \left[ \frac{\overline{T}_i-\overline{T}_i^{lev}}{\sigma_i^T}\right]^2
\end{equation}
where $\overline{T}_i$ and $\overline{T}_i^{lev}$ are, respectively, the model
and observed potential temperature at each
grid point $i$. The differences are weighted by an {\it a priori} uncertainty
$\sigma_i^T$ on observations (as provided by \cite{lev:94a}). The error
$\sigma_i^T$ is only a function of depth and varies from 0.5 at the surface to
0.05~K at the bottom of the ocean, mainly reflecting the decreasing
temperature variance with depth (Fig. \ref{Error}a). A value of $J_1$ of
order 1 means that the model is, on average, within observational
uncertainties.

The cost function also places constraints on the adjustment to insure it is
"reasonable", i.e. of order of the uncertainties on the observed surface heat
flux: 
\begin{equation}
J_2 = \frac{1}{N} \sum_{i=1}^N \left[\frac{Q_\mathrm{netm}}{\sigma^Q_i} \right]^2
\end{equation}
where $\sigma^Q_i$ are the {\it a priori} errors on the observed heat flux as
estimated by Stammer et al. (2002) from 30\% of local root-mean-square
variability of the NCEP forcing field (Fig \ref{Error}b).

The total cost function is defined as $J=\lambda_1 J_1+ \lambda_2 J_2$ where
$\lambda_1$ and $\lambda_2$ are weights controlling the relative contribution
of the two components. The adjoint model then yields the sensitivities
$\partial J/\partial Q_\mathrm{netm}$ of $J$ relative to the 2D fields
$Q_\mathrm{netm}$. Using a line-searching algorithm (\cite{gil-lem:89}),
$Q_\mathrm{netm}$ is adjusted then in the sense to
reduce $J$ --- the procedure is repeated until convergence.

%The configuration is identical 
%to the ``Global Ocean circulation'' tutorial where more details can be found.

Fig. \ref{Results} shows the results of such an optimization. The
model is started from rest and from January-mean temperature and salinity
initial conditions taken from the Levitus dataset. The experiment is run a year
and the averaged temperature over the whole run (i.e. annual mean) is used
in the cost function (\ref{cost_temp}) to evaluate the model\footnote{Because of
the daily automatic testing, the experiment as found in the repository is set-up
with a very small number of time-steps. To reproduce the results shown here, one
needs to set nTimeSteps = 360 and lastinterval=31104000 (both corresponding to a
year, see section \ref{sectRun} for further details).}. Only the
top 2 levels are used. The first guess $Q_\mathrm{netm}$ is chosen to be
zero. The weights $\lambda_1$ and $\lambda_2$ are set to 1 and 2, respectively.
The total cost function converges after 15 iterations, decreasing from 6.1 to
2.7 (the temperature contribution decreases from 6.1 to 1.8 while the heat
flux one increases from 0 to 0.42). The right panels of Fig. (\ref{Results})
illustrate the evolution of the temperature error at the surface from
iteration 0 to iteration 15. Unsurprisingly, the largest errors at iteration 0
(up to 6$^\circ$C, top left panels) are found in the Western boundary
currents. After optimization, the departure of the model temperature from
observations is reduced to 1$^\circ$C or less almost everywhere except in the
Pacific Equatorial Cold Tongue. Comparison of the initial temperature
error (top, right) and heat flux adjustment (bottom, left) shows that the
system basically increased the heat flux out of the ocean where temperatures
were too warm and vice-versa. Obviously, heat flux uncertainties are not the
sole responsible for temperature errors and the heat flux adjustment partly
compensates the poor representation of narrow currents (Western boundary
currents, Equatorial currents) at $4\times4^\circ$ resolution. This is
allowed by the large {\it a priori} error on the heat flux (Fig. \ref{Error}).
The Pacific Cold Tongue is a counter example: there, heat fluxes uncertainties
are fairly small (about 20~W.m$^2$), and a large temperature errors 
remains after optimization.

In the following, section 2 describes in details the implementation of the 
control variable $Q_\mathrm{netm}$, the cost function $J$ and the I/O required
for the communication between the model and the line-search. Instructions to
compile the MITgcm and its adjoint and the line-search algorithm are given in
section 3. The method used to run the experiment is described in section 4.

\begin{figure} [tpb]
\begin{center}
\includegraphics[width=\textwidth,height=.3\textheight]{s_examples/global_oce_optim/Error.eps}
\caption{{\it A priori} errors on potential temperature (left, in $^\circ$C) and
surface heat flux (right, in W~m$^{-2}$) used to compute the cost
terms $J_1$ and $J_2$, respectively.} 
\label{Error}
\end{center}
\end{figure}

\begin{figure} [tpb]
\begin{center}
\includegraphics[width=\textwidth,height=.3\textheight]{s_examples/global_oce_optim/Tutorial_fig.eps}
\caption{Initial annual mean surface heat flux (top right in W.m$^{-2}$) and 
adjustment obtained at iteration 15 (bottom right). Averaged difference
between model and observed potential temperatures at the surface (in $^\circ$C) 
before optimization (iteration 0, top right) and after optimization
(iteration 15, bottom right). Contour intervals for heat flux and temperature
are 25~W.m$^{-2}$ and 1$^\circ$C, respectively. A positive flux is out of the
ocean.}
\label{Results}
\end{center}
\end{figure}

\subsection{Implementation of the control variable and the cost function}

One of the goal of this tutorial is to illustrate how to implement a new
control variable. Most of this is fairly generic and is done in the ctrl
and cost packages found in the pkg/ directory. The modifications can be
tracked by the CPP option ALLOW\_HFLUXM\_CONTROL or the comment
cHFLUXM\_CONTROL. The more specific modifications required for the experiment
are found in verification/tutorial\_global\_oce\_optim/code\_ad. Here follows
a brief description of the implementation.

\subsubsection{The control variable}

The adjustment $Q_\mathrm{netm}$ is activated by setting
ALLOW\_HFLUXM\_CONTROL to "define" in ECCO\_OPTIONS.h.

It is first implemented as a ``normal'' forcing variable. It is defined in
FFIELDS.h, initialized to zero in ini\_forcing.F, and then used in
external\_forcing\_surf.F. $Q_\mathrm{netm}$ is made a control variable in
the ctrl package by modifying the following subroutines:

\begin{itemize}
\item ctrl\_init.F where $Q_\mathrm{netm}$ is defined as the control variable
number 24,

\item ctrl\_pack.F which writes, at the end of each iteration, the sensitivity
of the cost function $\partial J/\partial Q_\mathrm{netm}$ in to a file to be
used by the line-search algorithm,

\item ctrl\_unpack.F which reads, at the start of each iteration, the updated
adjustment as provided by the line-search algorithm,

\item ctrl\_map\_forcing.F in which the updated adjustment is added to the
first guess $Q_\mathrm{netm}$.
\end{itemize}

Note also some minor changes in ctrl.h, ctrl\_readparams.F, and ctrl\_dummy.h
(xx\_hfluxm\_file, fname\_hfluxm, xx\_hfluxm\_dummy).

\subsubsection{Cost functions}

The cost functions are implemented using the {\it cost} package.

\begin{itemize}

\item The temperature cost function $J_1$ which measures the drift of the mean
model temperature from the Levitus climatology is implemented in cost\_temp.F.
It is activated by ALLOW\_COST\_TEMP in ECCO\_OPTIONS.h. It requires the mean
temperature of the model which is obtained by accumulating the temperature in
cost\_tile.F (called at each time step).
The value of the cost function is stored in {\it objf\_temp} and its weight
$\lambda_1$ in {\it mult\_temp}.

\item The heat flux cost function, penalizing the departure of the surface
heat flux from observations is implemented in cost\_hflux.F, and activated by
the key ALLOW\_COST\_HFLUXM in ECCO\_OPTIONS.h. The value of the cost
function is stored in {\it objf\_hfluxm} and its weight $\lambda_2$ in
{\it mult\_hfluxm}.

\item The subroutine cost\_final.F calls the cost\_functions subroutines
and make the (weighted) sum of the various contributions.

\item The various weights used in the cost functions are read in
cost\_weights.F. The weight of the cost functions are read in
cost\_readparams.F from the input file data.cost.    

\end{itemize}


\subsection{Code Configuration}
%\label{www:tutorials}
\label{sec:eg_globest_code_config}

The model configuration for this experiment resides under the directory
{\it verification/tutorial\_global\_oce\_optim/}.  The experiment files in
code\_ad/ and input\_ad/ contain the code customizations and parameter
settings. Most of them are identical to those used in the Global Ocean
( experiment {\it  tutorial\_global\_oce\_latlon}). Below, we describe some of
the customizations required for this experiment.

\subsubsection{Compilation-time customizations in {\it code\_ad/}}

In ECCO\_CPPOPTIONS.h:

\begin{itemize}
\item define ALLOW\_ECCO\_OPTIMIZATION 

\item define ALLOW\_COST, ALLOW\_COST\_TEMP, and ALLOW\_COST\_HFLUXM

\item define ALLOW\_HFLUXM\_CONTROL
\end{itemize}

\subsubsection{Running-time customizations in {\it input\_ad/}}\label{sectRun}

\begin{itemize}

\item {\it data}: note the smaller {\it cg2dTargetResidual} than in the
forward-only experiment,

\item {\it data.optim} specifies the iteration number,

\item {\it data.ctrl} is used, in particular, to specify the
name of the sensitivity and adjustment files associated to a control
variable,

\item {\it data.cost}: parameters of the cost functions, in particular
{\it lastinterval} specifies the length of time-averaging for the model
temperature to be used in the cost function (\ref{cost_temp}),

\item {\it data.pkg}: note that the Gradient Check package is turned on by
default (useGrdchk=.TRUE.),

\item {\it Err\_hflux.bin} and {\it Err\_levitus\_15layer.bin} are the
files containing the heat flux and potential temperature uncertainties,
respectively.

\end{itemize}

\subsection{Compiling} 

The optimization experiment requires two executables: 1) the 
MITgcm and its adjoint ({\it mitgcmuv\_ad}) and 2) the line-search
algorithm ({\it optim.x}).

\subsubsection{Compilation of MITgcm and its adjoint: {\it mitcgmuv\_ad}}

Before compiling, first note that, in the directory code\_ad/, two files
must be updated:
\begin{itemize}
\item code\_ad\_diff.list which lists new subroutines to be compiled
by the TAF software (cost\_temp.F and cost\_hflux.F here),

\item  the adjoint\_hfluxm files which provides a list of the control variables
and the name of cost function to the TAF software.

\end{itemize}

Then, in the directory build\_ad/, type:
\begin{verbatim}
% ../../../tools/genmake2 -mods=../code\_ad -adof=../code\_ad/adjoint\_hfluxm
% make depend
% make adall
\end{verbatim}
to generate the MITgcm executable {\it mitgcmuv\_ad}.

\subsubsection{Compilation of the line-search algorithm: {\it optim.x}}

This is done from the directories lsopt/ and optim/ (under MITgcm/). In lsopt/,
unzip the blash1 library adapted to your platform, and change the Makefile
accordingly. Compile with:
\begin{verbatim}
% make all
\end{verbatim}
(more details in lsopt\_doc.txt)

In optim/,  the path of the directory where {\it mitgcm\_ad} was compiled
must be specified in the Makefile in the variable INCLUDEDIRS. The file name
of the control variable (xx\_hfluxm\_file here) must be added to the name list
read by optim\_num.F. Then use
\begin{verbatim}
% make depend
\end{verbatim} 
and
\begin{verbatim}
% make
\end{verbatim}
to generate the line-search executable {\it optim.x}.

\subsection{Running the estimation}

Copy the {\it mitgcmuv\_ad} executable to input\_ad/ and {\it optim.x} to the
subdirectory input\_ad/OPTIM/. Move into input\_ad/. The first iteration is
somewhat particular and is best done "by hand" while the following iterations
can be run automatically (see below). Check that the iteration number is set
to 0 in data.optim and run the MITgcm:
\begin{verbatim}
% ./mitgcmuv_ad
\end{verbatim}

The output files adxx\_hfluxm.0000000000.* and xx\_hfluxm.0000000000.* contain
the sensitivity of the cost function to $Q_\mathrm{netm}$ and the adjustment
to $Q_\mathrm{netm}$ (zero at the first iteration), respectively. Two other
files called costhflux\_tut\_MITgcm.opt0000 and ctrlhflux\_tut\_MITgcm.opt0000
are also generated. They essentially contain the same information as the
adxx\_.hfluxm* and xx\_hfluxm* files, but in a compressed format. These two files
are the only ones involved in the communication between the adjoint model
{\it mitgcmuv\_ad} and the line-search algorithm {\it optim.x}. Only at the first
iteration, are they both generated by {\it mitgcmuv\_ad}. Subsenquently,
costhflux\_tut\_MITgcm.opt$n$ is an output of the adjoint model at
iteration $n$ and an input of the line-search. The latter returns an updated
adjustment in ctrlhflux\_tut\_MITgcm.opt$n+1$ to be used as an input of the
adjoint model at iteration n+1. 

At the first iteration, move costhflux\_tut\_MITgcm.opt0000 and
ctrlhflux\_tut\_MITgcm.opt0000 to OPTIM/, move into this directory and link data.optim 
and data.ctrl locally:
\begin{verbatim}
% cd OPTIM/
% ln -s ../data.optim .
% ln -s ../data.ctrl .
\end{verbatim}
The target cost function {\it fmin} needs to be specified too: as a rule of thumb,
it should be about 0.95-0.90 times the value of the cost function at
the first iteration. This value is only used at the first iteration and does
not need to be updated afterward. However, it implicitly specifies the
``pace'' at which the cost function is going down (if you are lucky and it does
indeed diminish!). More in the ECCO section maybe ?

Once this is done, run the line-search algorithm:
\begin{verbatim}
% ./optim.x
\end{verbatim}
which computes the updated adjustment for iteration 1, ctrlhflux\_tut\_MITgcm.opt0001.

The following iterations can be executed automatically using the shell
script {\it cycsh} found in input\_ad/. This script will take care of changing
the iteration numbers in the data.optim, launch the adjoint model, clean and
store the outputs, move the costhflux* and ctrlhflux* files, and run the
line-search algorithm. Edit {\it cycsh} to specify the prefix of the
directories used to store the outputs and the maximum number of iteration.

1	jmc	1.17	% $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	dfer	1.1	% $Name: $
3
4	dfer	1.5	\section[Global Ocean State Estimation Example]{Global Ocean State Estimation at $4^\circ$ Resolution}
5	jmc	1.17	%\label{www:tutorials}
6			\label{sec:eg-global_state_estimate}
7	dfer	1.1	\begin{rawhtml}
8			<!-- CMIREDIR:eg-global_state_estimate: -->
9			\end{rawhtml}
10	jmc	1.9	\begin{center}
11			(in directory: {\it verification/tutorial\_global\_oce\_optim/})
12			\end{center}
13	dfer	1.1
14	dfer	1.3	\subsection{Overview}
15
16	dfer	1.12	This experiment illustrates the optimization capacity of the MITgcm: here,
17			a high level description.
18	dfer	1.1
19	dfer	1.12	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	dfer	1.1
27	dfer	1.10	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	dfer	1.12	method (see, for example, \cite{stam-etal:02} and \cite{fer-eta:05}).
31	dfer	1.1
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	dfer	1.12	\begin{equation}\label{cost_temp}
36	dfer	1.1	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	dfer	1.12	where $\overline{T}_i$ and $\overline{T}_i^{lev}$ are, respectively, the model
39			and observed potential temperature at each
40	dfer	1.10	grid point $i$. The differences are weighted by an {\it a priori} uncertainty
41	dfer	1.12	$\sigma_i^T$ on observations (as provided by \cite{lev:94a}). The error
42	dfer	1.10	$\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	dfer	1.1
48	dfer	1.10	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	dfer	1.1	flux:
51			\begin{equation}
52	dfer	1.10	J_2 = \frac{1}{N} \sum_{i=1}^N \left[\frac{Q_\mathrm{netm}}{\sigma^Q_i} \right]^2
53	dfer	1.1	\end{equation}
54	dfer	1.10	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	dfer	1.12	variability of the NCEP forcing field (Fig \ref{Error}b).
57	dfer	1.1
58	dfer	1.10	The total cost function is defined as $J=\lambda_1 J_1+ \lambda_2 J_2$ where
59	dfer	1.1	$\lambda_1$ and $\lambda_2$ are weights controlling the relative contribution
60	dfer	1.10	of the two components. The adjoint model then yields the sensitivities
61	dfer	1.1	$\partial J/\partial Q_\mathrm{netm}$ of $J$ relative to the 2D fields
62	dfer	1.12	$Q_\mathrm{netm}$. Using a line-searching algorithm (\cite{gil-lem:89}),
63			$Q_\mathrm{netm}$ is adjusted then in the sense to
64	dfer	1.10	reduce $J$ --- the procedure is repeated until convergence.
65
66	dfer	1.12	%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	dfer	1.15	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	dfer	1.12	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	dfer	1.10
106			\begin{figure} [tpb]
107			\begin{center}
108	jmc	1.16	\includegraphics[width=\textwidth,height=.3\textheight]{s_examples/global_oce_optim/Error.eps}
109	dfer	1.12	\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	dfer	1.10	\label{Error}
113			\end{center}
114			\end{figure}
115	dfer	1.1
116	dfer	1.12	\begin{figure} [tpb]
117			\begin{center}
118	jmc	1.16	\includegraphics[width=\textwidth,height=.3\textheight]{s_examples/global_oce_optim/Tutorial_fig.eps}
119	dfer	1.12	\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	dfer	1.1
130	dfer	1.10	\subsection{Implementation of the control variable and the cost function}
131	dfer	1.1
132	dfer	1.12	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	dfer	1.1
140	dfer	1.3	\subsubsection{The control variable}
141
142	dfer	1.10	The adjustment $Q_\mathrm{netm}$ is activated by setting
143			ALLOW\_HFLUXM\_CONTROL to "define" in ECCO\_OPTIONS.h.
144	dfer	1.1
145	dfer	1.10	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	dfer	1.1
150	dfer	1.3	\begin{itemize}
151	dfer	1.10	\item ctrl\_init.F where $Q_\mathrm{netm}$ is defined as the control variable
152			number 24,
153	dfer	1.1
154	dfer	1.10	\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	dfer	1.1
158	dfer	1.10	\item ctrl\_unpack.F which reads, at the start of each iteration, the updated
159	dfer	1.12	adjustment as provided by the line-search algorithm,
160	dfer	1.1
161	dfer	1.12	\item ctrl\_map\_forcing.F in which the updated adjustment is added to the
162	dfer	1.10	first guess $Q_\mathrm{netm}$.
163	dfer	1.3	\end{itemize}
164	dfer	1.1
165	dfer	1.10	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	dfer	1.1
168	dfer	1.3	\subsubsection{Cost functions}
169
170	dfer	1.10	The cost functions are implemented using the {\it cost} package.
171	dfer	1.1
172	dfer	1.3	\begin{itemize}
173	dfer	1.1
174	dfer	1.10	\item The temperature cost function $J_1$ which measures the drift of the mean
175	dfer	1.12	model temperature from the Levitus climatology is implemented in cost\_temp.F.
176	dfer	1.10	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	dfer	1.12	the key ALLOW\_COST\_HFLUXM in ECCO\_OPTIONS.h. The value of the cost
185	dfer	1.10	function is stored in {\it objf\_hfluxm} and its weight $\lambda_2$ in
186			{\it mult\_hfluxm}.
187	dfer	1.1
188	dfer	1.3	\item The subroutine cost\_final.F calls the cost\_functions subroutines
189	dfer	1.10	and make the (weighted) sum of the various contributions.
190	dfer	1.1
191	dfer	1.10	\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	dfer	1.1
195	dfer	1.3	\end{itemize}
196
197	dfer	1.5
198			\subsection{Code Configuration}
199	jmc	1.17	%\label{www:tutorials}
200			\label{sec:eg_globest_code_config}
201	dfer	1.5
202	dfer	1.10	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	jmc	1.13	( experiment {\it tutorial\_global\_oce\_latlon}). Below, we describe some of
207	dfer	1.12	the customizations required for this experiment.
208	dfer	1.10
209			\subsubsection{Compilation-time customizations in {\it code\_ad/}}
210
211	dfer	1.12	In ECCO\_CPPOPTIONS.h:
212	dfer	1.10
213			\begin{itemize}
214	dfer	1.12	\item define ALLOW\_ECCO\_OPTIMIZATION
215	dfer	1.10
216	dfer	1.12	\item define ALLOW\_COST, ALLOW\_COST\_TEMP, and ALLOW\_COST\_HFLUXM
217	dfer	1.10
218	dfer	1.12	\item define ALLOW\_HFLUXM\_CONTROL
219	dfer	1.10	\end{itemize}
220
221	dfer	1.15	\subsubsection{Running-time customizations in {\it input\_ad/}}\label{sectRun}
222	dfer	1.10
223			\begin{itemize}
224
225	dfer	1.12	\item {\it data}: note the smaller {\it cg2dTargetResidual} than in the
226			forward-only experiment,
227	dfer	1.10
228	dfer	1.12	\item {\it data.optim} specifies the iteration number,
229	dfer	1.10
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	dfer	1.12	variable,
233	dfer	1.10
234	dfer	1.12	\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	dfer	1.10
238	dfer	1.12	\item {\it data.pkg}: note that the Gradient Check package is turned on by
239			default (useGrdchk=.TRUE.),
240	dfer	1.10
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	dfer	1.5
245	dfer	1.10	\end{itemize}
246	dfer	1.5
247	dfer	1.3	\subsection{Compiling}
248	dfer	1.1
249			The optimization experiment requires two executables: 1) the
250	dfer	1.10	MITgcm and its adjoint ({\it mitgcmuv\_ad}) and 2) the line-search
251			algorithm ({\it optim.x}).
252	dfer	1.1
253	edhill	1.6	\subsubsection{Compilation of MITgcm and its adjoint: {\it mitcgmuv\_ad}}
254	dfer	1.1
255			Before compiling, first note that, in the directory code\_ad/, two files
256			must be updated:
257	dfer	1.5	\begin{itemize}
258			\item code\_ad\_diff.list which lists new subroutines to be compiled
259	dfer	1.10	by the TAF software (cost\_temp.F and cost\_hflux.F here),
260	dfer	1.1
261	dfer	1.5	\item the adjoint\_hfluxm files which provides a list of the control variables
262	dfer	1.10	and the name of cost function to the TAF software.
263	dfer	1.5
264			\end{itemize}
265	dfer	1.1
266			Then, in the directory build\_ad/, type:
267	dfer	1.5	\begin{verbatim}
268			% ../../../tools/genmake2 -mods=../code\_ad -adof=../code\_ad/adjoint\_hfluxm
269			% make depend
270			% make adall
271			\end{verbatim}
272	dfer	1.10	to generate the MITgcm executable {\it mitgcmuv\_ad}.
273	dfer	1.1
274	edhill	1.6	\subsubsection{Compilation of the line-search algorithm: {\it optim.x}}
275	dfer	1.1
276	dfer	1.10	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	dfer	1.5	\begin{verbatim}
280			% make all
281			\end{verbatim}
282			(more details in lsopt\_doc.txt)
283	dfer	1.1
284	dfer	1.10	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	dfer	1.5	\begin{verbatim}
289			% make depend
290			\end{verbatim}
291			and
292			\begin{verbatim}
293			% make
294			\end{verbatim}
295	edhill	1.6	to generate the line-search executable {\it optim.x}.
296	dfer	1.1
297	dfer	1.3	\subsection{Running the estimation}
298	dfer	1.1
299	dfer	1.10	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	dfer	1.3	\begin{verbatim}
305			% ./mitgcmuv_ad
306			\end{verbatim}
307	dfer	1.1
308	dfer	1.5	The output files adxx\_hfluxm.0000000000.* and xx\_hfluxm.0000000000.* contain
309	dfer	1.12	the sensitivity of the cost function to $Q_\mathrm{netm}$ and the adjustment
310	dfer	1.10	to $Q_\mathrm{netm}$ (zero at the first iteration), respectively. Two other
311	dfer	1.12	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	dfer	1.10	iteration $n$ and an input of the line-search. The latter returns an updated
319	dfer	1.12	adjustment in ctrlhflux\_tut\_MITgcm.opt$n+1$ to be used as an input of the
320			adjoint model at iteration n+1.
321	dfer	1.10
322	dfer	1.12	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	dfer	1.10	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	dfer	1.4
337			Once this is done, run the line-search algorithm:
338			\begin{verbatim}
339			% ./optim.x
340			\end{verbatim}
341	dfer	1.12	which computes the updated adjustment for iteration 1, ctrlhflux\_tut\_MITgcm.opt0001.
342	dfer	1.4
343	edhill	1.6	The following iterations can be executed automatically using the shell
344	dfer	1.10	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	dfer	1.12	store the outputs, move the costhflux* and ctrlhflux* files, and run the
347	dfer	1.10	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	dfer	1.12