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.3 2005/08/02 19:54:53 dfer 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}

\subsection{Overview}

Mean surface heat flux as a  control variable : Qnetm

This experiment illustrates the optimization (or data-assimilation) capacity
of the MITgcm. Using an ocean configuration with realistic geography and bathymetry on a
4x4 spherical polar grid, we estimate a time-independent surface heat flux correction
Qnetm that brings the model climatology into consistency with observations (Levitus
climatology).

This correction Qnetm (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 Qnetm 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 Qnetm 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.

Qnetm is made a control variable in the ctrl package by modifying the following subroutines:

\begin{itemize}
\item ctrl\_init.F where Qnetm 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 Qnetm.
\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{Compiling} 

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

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

Before compiling, first note that, in the directory code\_ad/, two files
must be updated:
- the code\_ad\_diff.list file which lists new subroutines which are to be compiled
by the TAF software (cost\_temp.f and cost\_hflux.f here),

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

Then, in the directory build\_ad/, type:
../../../tools/genmake2 -mods=../code\_ad -adof=../code\_ad/adjoint\_hfluxm
make depend
make adall
to generate the MITgcm excutable 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 blas1 library you need, and change accordingly the
library in the Makefile. Compile with make all

(more details in lsopt\_doc.txt)
In optim/,  the path of the directory where 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 \it{make depend} and \it{make} to generate the line-search executable \it{optim.x}.

\subsection{Running the estimation}

Copy the mitgcmuv\_ad executble to input\_ad and 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,
respectively, the sensitivity of the cost function to Qnetm and the perturbations
to Qnetm (zero at the first iteration). 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 produces the updated perturbation for iterations 1 ecco\_ctrl\_1.

The following iterations can be executed automatically using the shell script \it{cycle}
found in input\_ad/.

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