1 |
% $Header: /u/gcmpack/manual/part3/case_studies/climatalogical_ogcm/climatalogical_ogcm.tex,v 1.12 2004/10/16 03:40:13 edhill Exp $ |
2 |
% $Name: $ |
3 |
|
4 |
\section[Global Ocean State Estimation Example]{Global Ocean State Estimation at 4$^\circ$ Resolution} |
5 |
\label{www:tutorials} |
6 |
\label{sect:eg-global_state_estimate} |
7 |
\begin{rawhtml} |
8 |
<!-- CMIREDIR:eg-global_state_estimate: --> |
9 |
\end{rawhtml} |
10 |
|
11 |
1) Overview: Mean surface heat flux as a control variable : Qnetm |
12 |
|
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. |