s_examples/global_oce_optim/global_oce_estimation.tex

% $Header: /u/gcmpack/manual/part3/case_studies/global_oce_estimation/global_oce_estimation.tex,v 1.8 2006/06/28 19:00:53 molod Exp $
% $Name:  $

\section[Global Ocean State Estimation Example]{Global Ocean State Estimation at $4^\circ$ Resolution}
\label{www:tutorials}
\label{sect: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}

Mean surface heat flux as a  control variable : $Q_\mathrm{netm}$

This experiment illustrates the optimization (or data-assimilation) 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 correction
$Q_\mathrm{netm}$ that brings the model climatology into consistency with observations (Levitus
climatology). The files for this experiment can be found in the verification directory under
tutorial\_global\_oce\_optim.

This correction $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, 
Stammer et al. (2002) and Ferreira et al. (2005)).

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}
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$ is the averaged model temperature and $\overline{T}_i^{lev}$
the annual mean observed temperature at each grid point $i$. The differences
are weighted by an a priori uncertainty $\sigma_i^T$ on observations (Levitus
and Boyer 1994). 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. 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 correction to insure it is
"reasonnable", 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^x_i} \right]^2
\end{equation}
where $\sigma^x_i$ are the a priori errors (2d field from ECCO ..... Fig ?).

The total cost function is obtained 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 mcomponents. The adjoint model then provides the sensitivities
$\partial J/\partial Q_\mathrm{netm}$ of $J$ relative to the 2D fields
$Q_\mathrm{netm}$. Using a line-searching algorithm (Gilbert and Lemar\'{e}chal 1989),
$Q_\mathrm{netm}$ is adjusted in the sense to reduce $J$ --- the procedure is 
repeated until convergence.

In the following example, the configuration is identical to the "Global ocean circulation"
tutorial where more details can be found. In each iteration, the model is started from
rest with temperature and salinity initial conditions taken from Levitus dataset and run
for a year. The first guess $Q_\mathrm{netm}$ is chosen to be zero.

The experiment employs two executables: one for the MITgcm and its adjoints and 
one for the line-search algorithmi (offline optimization). The implementation of
the control variable $Q_\mathrm{netm}$, the cost function $J$ and the I/O required
for the commmunication betwwen the model and the line-search are described in details 
in section 2. The compilation of the two executables is given in section 3.
A method to run the experiment is described in section 4.

Gilbert, J. C., and C. Lemar\'echal, 1989: Some numerical experiments with
variable-storage quasi-Newton algorithms. \textit{Math. Programm.,}
\textbf{45,} 407-435.

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

All subroutines that require modifications are found in verifications/Optim/code\_ad

\subsubsection{The control variable}

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

It is first implemented as a 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 iteration, the sensitivity of the cost function
$\partial J/\partial Q_\mathrm{netm}$ into a file to be used by the lins-search algorithm,

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

\item ctrl\_map\_forcing.F where the updated perturbation 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 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 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 objf\_temp and its weight $\lambda_1$
in 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 activated
ALLOW\_COST\_HFLUXM in ECCO\_OPTIONS.h. The value of the cost function is stored in
objf\_hfluxm and its weight $\lambda_2$ in mult\_hfluxm,

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

\item The weights used in the cost functions are read is cost\_weights.F.
The weigth 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_fourl_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 customisations and parameter settings
for this experiment. Most of them are identical to those used in
the Global ocean experiment. Below we describe the customisations to
these files associated with this experiment.

\subsubsection{File {\it input/data}}


\subsection{Compiling} 

The optimization experiment requires two executables: 1) the 
MITgcm and its adjoints (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 sotware.

\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 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 you need, and change accordingly the
library in the Makefile. 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 compliled must be specified
in the Makefile in the variable INCLUDEDIRS. The file name of the controle variable
(xx\_hfluxm\_file here) must be added to the namelist 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 has be done "by hand". 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 perturbations
to $Q_\mathrm{netm}$ (zero at the first iteration), respectively. Two other files
called ecco\_cost.. and ecco\_ctrl are also generated. They essentially contains
the same information as the adxx\_* and xx\_* files, but in a compressed format.
These two file are the only ones involved in the communication between the adjoint
model {\it mitgcmuv\_ad} and the line-search algorithm {\it optim.x}. The ecco\_cost*n
is an ouput of the adjoint model at iteration $n$ and an input of the line-search. The
latter returns an updated perturbation in ecco\_ctrl*n+1 to be used as an input of
the adjoint model at iteration n+1. 

At the first iteration, move these two files ecco\_cost and ecco\_ctrl 
to OPTIM/, open data.optim and check the iteration number is set to 0.
The target cost function {\it fmin} needs to be specified: a rule of thumb 
suggest it should be set to 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 afterwards. 
However, it implicitly specifies the "pace" at which the cost function is
going down (if you are lucky and it goes down). 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 perturbations for iteration 1 ecco\_ctrl\_1.

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 ouputs, move the
ecco\_cost* and ecco\_ctrl* files, and run the line-search algotrithm.
Edit {\it cycsh} to specify the prefix of the directories used to store the ouputs and
the maximum number of iteration.

1	jmc	1.9	% $Header: /u/gcmpack/manual/part3/case_studies/global_oce_estimation/global_oce_estimation.tex,v 1.8 2006/06/28 19:00:53 molod 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	dfer	1.1	\label{www:tutorials}
6			\label{sect:eg-global_state_estimate}
7			\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.5	Mean surface heat flux as a control variable : $Q_\mathrm{netm}$
17	dfer	1.1
18			This experiment illustrates the optimization (or data-assimilation) capacity
19			of the MITgcm. Using an ocean configuration with realistic geography and bathymetry on a
20	dfer	1.5	$4\times4^circ$ spherical polar grid, we estimate a time-independent surface heat flux correction
21			$Q_\mathrm{netm}$ that brings the model climatology into consistency with observations (Levitus
22	molod	1.7	climatology). The files for this experiment can be found in the verification directory under
23			tutorial\_global\_oce\_optim.
24	dfer	1.1
25	dfer	1.5	This correction $Q_\mathrm{netm}$ (a 2D field only function of longitude and latitude) is
26	dfer	1.1	the control variable of an optimization problem. It is inferred by an iterative
27			procedure using an `adjoint technique' and a least-squares method (see, for example,
28			Stammer et al. (2002) and Ferreira et al. (2005)).
29
30			The ocean model is run forward in time and the quality of the solution is
31			determined by a cost function, $J_1$, a measure of the departure of the model
32			climatology from observations:
33			\begin{equation}
34			J_1=\frac{1}{N}\sum_{i=1}^N \left[ \frac{\overline{T}_i-\overline{T}_i^{lev}}{\sigma_i^T}\right]^2
35			\end{equation}
36			where $\overline{T}_i$ is the averaged model temperature and $\overline{T}_i^{lev}$
37			the annual mean observed temperature at each grid point $i$. The differences
38			are weighted by an a priori uncertainty $\sigma_i^T$ on observations (Levitus
39			and Boyer 1994). The error $\sigma_i^T$ is only a function of depth and varies
40			from 0.5 at the surface to 0.05~K at the bottom of the ocean, mainly reflecting
41			the decreasing temperature variance with depth. A value of $J_1$ of order 1 means
42			that the model is, on average, within observational uncertainties.
43
44			The cost function also places constraints on the correction to insure it is
45			"reasonnable", i.e. of order of the uncertainties on the observed surface heat
46			flux:
47			\begin{equation}
48			J_2 = \frac{1}{N} \sum_{i=1}^N \left[\frac{Q_\mathrm{netm}}{\sigma^x_i} \right]^2
49			\end{equation}
50			where $\sigma^x_i$ are the a priori errors (2d field from ECCO ..... Fig ?).
51
52			The total cost function is obtained as $J=\lambda_1 J_1+ \lambda_2 J_2$ where
53			$\lambda_1$ and $\lambda_2$ are weights controlling the relative contribution
54			of the two mcomponents. The adjoint model then provides the sensitivities
55			$\partial J/\partial Q_\mathrm{netm}$ of $J$ relative to the 2D fields
56			$Q_\mathrm{netm}$. Using a line-searching algorithm (Gilbert and Lemar\'{e}chal 1989),
57			$Q_\mathrm{netm}$ is adjusted in the sense to reduce $J$ --- the procedure is
58			repeated until convergence.
59
60			In the following example, the configuration is identical to the "Global ocean circulation"
61			tutorial where more details can be found. In each iteration, the model is started from
62			rest with temperature and salinity initial conditions taken from Levitus dataset and run
63	dfer	1.5	for a year. The first guess $Q_\mathrm{netm}$ is chosen to be zero.
64	dfer	1.1
65			The experiment employs two executables: one for the MITgcm and its adjoints and
66			one for the line-search algorithmi (offline optimization). The implementation of
67			the control variable $Q_\mathrm{netm}$, the cost function $J$ and the I/O required
68			for the commmunication betwwen the model and the line-search are described in details
69			in section 2. The compilation of the two executables is given in section 3.
70	dfer	1.3	A method to run the experiment is described in section 4.
71	dfer	1.1
72			Gilbert, J. C., and C. Lemar\'echal, 1989: Some numerical experiments with
73			variable-storage quasi-Newton algorithms. \textit{Math. Programm.,}
74			\textbf{45,} 407-435.
75
76	dfer	1.3	\subsection{Implemention of the control variable and the cost function}
77	dfer	1.1
78			All subroutines that require modifications are found in verifications/Optim/code\_ad
79
80	dfer	1.3	\subsubsection{The control variable}
81
82	dfer	1.5	The correction $Q_\mathrm{netm}$ is activated by setting ALLOW\_HFLUXM\_CONTROL to "define" in ECCO\_OPTIONS.h.
83	dfer	1.1
84			It is first implemented as a forcing variable. It is defined in FFIELDS.h,
85			initialized to zero in ini\_forcing.F, and then used in external\_forcing\_surf.F.
86
87	dfer	1.5	$Q_\mathrm{netm}$ is made a control variable in the ctrl package by modifying the following subroutines:
88	dfer	1.1
89	dfer	1.3	\begin{itemize}
90	dfer	1.5	\item ctrl\_init.F where $Q_\mathrm{netm}$ is defined as the control variable number 24,
91	dfer	1.1
92	dfer	1.3	\item ctrl\_pack.F which writes, at the end of iteration, the sensitivity of the cost function
93	dfer	1.1	$\partial J/\partial Q_\mathrm{netm}$ into a file to be used by the lins-search algorithm,
94
95	dfer	1.3	\item ctrl\_unpack.F which reads, at the start of each iteration, the updated perturbations as
96	dfer	1.1	provided by the line-search algorithm,
97
98	dfer	1.5	\item ctrl\_map\_forcing.F where the updated perturbation is added to the first guess $Q_\mathrm{netm}$.
99	dfer	1.3	\end{itemize}
100	dfer	1.1
101			Note also some minor changes in ctrl.h, ctrl\_readparams.F, and ctrl\_dummy.h (xx\_hfluxm\_file,
102			fname\_hfluxm, xx\_hfluxm\_dummy).
103
104	dfer	1.3	\subsubsection{Cost functions}
105
106			The cost functions are implemented using the cost package.
107	dfer	1.1
108	dfer	1.3	\begin{itemize}
109	dfer	1.1
110	dfer	1.3	\item The temperature cost function $J_1$ which measures the drift of the mean model
111	dfer	1.1	temperature from the Levitus climatology is implemented cost\_temp.F. It is
112			activated by ALLOW\_COST\_TEMP in ECCO\_OPTIONS.h. It requires the mean temperature of
113			the model which is obtained by accumulating the temperature in cost\_tile.F (called at
114			each time step).
115			The value of the cost function is stored in objf\_temp and its weight $\lambda_1$
116			in mult\_temp.
117
118	dfer	1.3	\item The heat flux cost function penalizing the departure of the surface heat flux from
119	dfer	1.1	observations is implemented in cost\_hflux.F, and activated by activated
120			ALLOW\_COST\_HFLUXM in ECCO\_OPTIONS.h. The value of the cost function is stored in
121			objf\_hfluxm and its weight $\lambda_2$ in mult\_hfluxm,
122
123	dfer	1.3	\item The subroutine cost\_final.F calls the cost\_functions subroutines
124	dfer	1.1	and make the (weighted) sum of the different contributions.
125
126	dfer	1.3	\item The weights used in the cost functions are read is cost\_weights.F.
127	dfer	1.1	The weigth of the cost functions are read in cost\_readparams.F from the input file
128			data.cost.
129
130	dfer	1.3	\end{itemize}
131
132	dfer	1.5
133			\subsection{Code Configuration}
134			\label{www:tutorials}
135			\label{SEC:eg_fourl_code_config}
136
137			The model configuration for this experiment resides under the
138	molod	1.8	directory {\it verification/tutorial\_global\_oce\_optim/}. The experiment files in code\_ad/
139	dfer	1.5	and input\_ad/ contain the code customisations and parameter settings
140			for this experiment. Most of them are identical to those used in
141			the Global ocean experiment. Below we describe the customisations to
142			these files associated with this experiment.
143
144			\subsubsection{File {\it input/data}}
145
146
147	dfer	1.3	\subsection{Compiling}
148	dfer	1.1
149			The optimization experiment requires two executables: 1) the
150	dfer	1.5	MITgcm and its adjoints (it{mitgcmuv\_ad}) and 2) the line-search
151	edhill	1.6	algorithm ({\it optim.x})
152	dfer	1.1
153	edhill	1.6	\subsubsection{Compilation of MITgcm and its adjoint: {\it mitcgmuv\_ad}}
154	dfer	1.1
155			Before compiling, first note that, in the directory code\_ad/, two files
156			must be updated:
157	dfer	1.5	\begin{itemize}
158			\item code\_ad\_diff.list which lists new subroutines to be compiled
159	dfer	1.1	by the TAF software (cost\_temp.f and cost\_hflux.f here),
160
161	dfer	1.5	\item the adjoint\_hfluxm files which provides a list of the control variables
162			and the name of cost function to the TAF sotware.
163
164			\end{itemize}
165	dfer	1.1
166			Then, in the directory build\_ad/, type:
167	dfer	1.5	\begin{verbatim}
168			% ../../../tools/genmake2 -mods=../code\_ad -adof=../code\_ad/adjoint\_hfluxm
169			% make depend
170			% make adall
171			\end{verbatim}
172			to generate the MITgcm executable mitgcmuv\_ad
173	dfer	1.1
174	edhill	1.6	\subsubsection{Compilation of the line-search algorithm: {\it optim.x}}
175	dfer	1.1
176	dfer	1.3	This is done from the directories lsopt/ and optim/ (under MITgcm/)
177	dfer	1.1
178	dfer	1.5	In lsopt/, unzip the blash1 library you need, and change accordingly the
179			library in the Makefile. Compile with:
180			\begin{verbatim}
181			% make all
182			\end{verbatim}
183			(more details in lsopt\_doc.txt)
184	dfer	1.1
185	edhill	1.6	In optim/, the path of the directory where {\it mitgcm\_ad} was compliled must be specified
186	dfer	1.1	in the Makefile in the variable INCLUDEDIRS. The file name of the controle variable
187			(xx\_hfluxm\_file here) must be added to the namelist read by optim\_num.F
188
189	dfer	1.5	Then use
190			\begin{verbatim}
191			% make depend
192			\end{verbatim}
193			and
194			\begin{verbatim}
195			% make
196			\end{verbatim}
197	edhill	1.6	to generate the line-search executable {\it optim.x}.
198	dfer	1.1
199	dfer	1.3	\subsection{Running the estimation}
200	dfer	1.1
201	edhill	1.6	Copy the {\it mitgcmuv\_ad} executable to input\_ad and {\it optim.x}
202	dfer	1.5	to the subdirectory input\_ad/OPTIM. Move into input\_ad/.
203	dfer	1.1
204	dfer	1.3	The first iteration has be done "by hand". Check that the iteration number is set
205			to 0 in data.optim and run the Mitgcm
206			\begin{verbatim}
207			% ./mitgcmuv_ad
208			\end{verbatim}
209	dfer	1.1
210	dfer	1.5	The output files adxx\_hfluxm.0000000000.* and xx\_hfluxm.0000000000.* contain
211			the sensitivity of the cost function to $Q_\mathrm{netm}$ and the perturbations
212			to $Q_\mathrm{netm}$ (zero at the first iteration), respectively. Two other files
213			called ecco\_cost.. and ecco\_ctrl are also generated. They essentially contains
214			the same information as the adxx\_* and xx\_* files, but in a compressed format.
215			These two file are the only ones involved in the communication between the adjoint
216	edhill	1.6	model {\it mitgcmuv\_ad} and the line-search algorithm {\it optim.x}. The ecco\_cost*n
217	dfer	1.5	is an ouput of the adjoint model at iteration $n$ and an input of the line-search. The
218			latter returns an updated perturbation in ecco\_ctrl*n+1 to be used as an input of
219			the adjoint model at iteration n+1.
220	dfer	1.4
221			At the first iteration, move these two files ecco\_cost and ecco\_ctrl
222			to OPTIM/, open data.optim and check the iteration number is set to 0.
223	edhill	1.6	The target cost function {\it fmin} needs to be specified: a rule of thumb
224	dfer	1.4	suggest it should be set to about 0.95-0.90 times the value of the cost
225			function at the first iteration. This value is only used at the first
226			iteration and does not need to be updated afterwards.
227			However, it implicitly specifies the "pace" at which the cost function is
228			going down (if you are lucky and it goes down). More in the ECCO section maybe ?
229
230			Once this is done, run the line-search algorithm:
231			\begin{verbatim}
232			% ./optim.x
233			\end{verbatim}
234	dfer	1.5	which computes the updated perturbations for iteration 1 ecco\_ctrl\_1.
235	dfer	1.4
236	edhill	1.6	The following iterations can be executed automatically using the shell
237			script {\it cycsh}
238	dfer	1.5	found in input\_ad/. This script will take care of changing the iteration numbers in the
239			data.optim, launch the adjoint model, clean and store the ouputs, move the
240			ecco\_cost* and ecco\_ctrl* files, and run the line-search algotrithm.
241	edhill	1.6	Edit {\it cycsh} to specify the prefix of the directories used to store the ouputs and
242	dfer	1.5	the maximum number of iteration.
243	dfer	1.1