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.1 2005/08/02 19:16:58 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}

1) 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.

2) Implemention of the control variable and the cost function.

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

2.1) 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:

- ctrl\_init.F where Qnetm is defined as the control variable number 24,

- 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,

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

- ctrl\_map\_forcing.F where the updated perturbation is added to the first guess Qnetm.

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

2.2) Cost functions

The cost functions are implemented in the cost package.

2.2.1) 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.

2.2.2) The 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,

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

2.2.4) 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.    

3) Compiling 

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

3.1) 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


3.2) Compilation of the line Search Algorithm (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 make depend and make to generate the line-search executable optim.x.

4) Running the optimization

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
./mitgcmuv\_ad

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, 
these are the only files involved in the communication between the model mitgcmuv\_ad and
the line search optim.x. the cost one contains the cost function and the cost function
sensitivity to the control variables. It is an ouptu of mitgcmuv\_ad. The ctrl one contain the correction to Qnetm.

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