JURASSIC
Functions
retrieval.c File Reference

JURASSIC retrieval processor. More...

#include "jurassic.h"

Go to the source code of this file.

Functions

void optimal_estimation (ret_t *ret, ctl_t *ctl, tbl_t *tbl, obs_t *obs_meas, obs_t *obs_i, atm_t *atm_apr, atm_t *atm_i)
 Carry out optimal estimation retrieval. More...
 
int main (int argc, char *argv[])
 

Detailed Description

JURASSIC retrieval processor.

Definition in file retrieval.c.

Function Documentation

◆ optimal_estimation()

void optimal_estimation ( ret_t *  ret,
ctl_t *  ctl,
tbl_t *  tbl,
obs_t *  obs_meas,
obs_t *  obs_i,
atm_t *  atm_apr,
atm_t *  atm_i 
)

Carry out optimal estimation retrieval.

Definition at line 107 of file retrieval.c.

114 {
115
116 static int ipa[N], iqa[N];
117
118 FILE *out;
119
120 char filename[2 * LEN];
121
122 double chisq, disq = 0, lmpar = 0.001;
123
124 int it = 0;
125
126 /* ------------------------------------------------------------
127 Initialize...
128 ------------------------------------------------------------ */
129
130 /* Get sizes... */
131 const size_t m = obs2y(ctl, obs_meas, NULL, NULL, NULL);
132 const size_t n = atm2x(ctl, atm_apr, NULL, iqa, ipa);
133 if (m == 0 || n == 0)
134 ERRMSG("Check problem definition!");
135
136 /* Write info... */
137 LOG(1, "Problem size: m= %d / n= %d "
138 "(alloc= %.4g MB / stat= %.4g MB)",
139 (int) m, (int) n,
140 (double) (3 * m * n + 4 * n * n + 8 * m +
141 8 * n) * sizeof(double) / 1024. / 1024.,
142 (double) (5 * sizeof(atm_t) + 3 * sizeof(obs_t)
143 + 2 * N * sizeof(int)) / 1024. / 1024.);
144
145 /* Allocate... */
146 gsl_matrix *a = gsl_matrix_alloc(n, n);
147 gsl_matrix *cov = gsl_matrix_alloc(n, n);
148 gsl_matrix *k_i = gsl_matrix_alloc(m, n);
149 gsl_matrix *s_a_inv = gsl_matrix_alloc(n, n);
150
151 gsl_vector *b = gsl_vector_alloc(n);
152 gsl_vector *dx = gsl_vector_alloc(n);
153 gsl_vector *dy = gsl_vector_alloc(m);
154 gsl_vector *sig_eps_inv = gsl_vector_alloc(m);
155 gsl_vector *sig_formod = gsl_vector_alloc(m);
156 gsl_vector *sig_noise = gsl_vector_alloc(m);
157 gsl_vector *x_a = gsl_vector_alloc(n);
158 gsl_vector *x_i = gsl_vector_alloc(n);
159 gsl_vector *x_step = gsl_vector_alloc(n);
160 gsl_vector *y_aux = gsl_vector_alloc(m);
161 gsl_vector *y_i = gsl_vector_alloc(m);
162 gsl_vector *y_m = gsl_vector_alloc(m);
163
164 /* Set initial state... */
165 copy_atm(ctl, atm_i, atm_apr, 0);
166 copy_obs(ctl, obs_i, obs_meas, 0);
167 formod(ctl, tbl, atm_i, obs_i);
168
169 /* Set state vectors and observation vectors... */
170 atm2x(ctl, atm_apr, x_a, NULL, NULL);
171 atm2x(ctl, atm_i, x_i, NULL, NULL);
172 obs2y(ctl, obs_meas, y_m, NULL, NULL);
173 obs2y(ctl, obs_i, y_i, NULL, NULL);
174
175 /* Set inverse a priori covariance S_a^-1... */
176 set_cov_apr(ret, ctl, atm_apr, iqa, ipa, s_a_inv);
177 write_matrix(ret->dir, "matrix_cov_apr.tab", ctl, s_a_inv,
178 atm_i, obs_i, "x", "x", "r");
179 matrix_invert(s_a_inv);
180
181 /* Get measurement errors... */
182 set_cov_meas(ret, ctl, obs_meas, sig_noise, sig_formod, sig_eps_inv);
183
184 /* Create cost function file... */
185 sprintf(filename, "%s/costs.tab", ret->dir);
186 if (!(out = fopen(filename, "w")))
187 ERRMSG("Cannot create cost function file!");
188
189 /* Write header... */
190 fprintf(out,
191 "# $1 = iteration number\n"
192 "# $2 = normalized cost function\n"
193 "# $3 = number of measurements\n"
194 "# $4 = number of state vector elements\n\n");
195
196 /* Determine dx = x_i - x_a and dy = y - F(x_i) ... */
197 gsl_vector_memcpy(dx, x_i);
198 gsl_vector_sub(dx, x_a);
199 gsl_vector_memcpy(dy, y_m);
200 gsl_vector_sub(dy, y_i);
201
202 /* Compute cost function... */
203 chisq = cost_function(dx, dy, s_a_inv, sig_eps_inv);
204
205 /* Write info... */
206 LOG(1, "it= %d / chi^2/m= %g", it, chisq);
207
208 /* Write to cost function file... */
209 fprintf(out, "%d %g %d %d\n", it, chisq, (int) m, (int) n);
210
211 /* Compute initial kernel... */
212 kernel(ctl, tbl, atm_i, obs_i, k_i);
213
214 /* ------------------------------------------------------------
215 Levenberg-Marquardt minimization...
216 ------------------------------------------------------------ */
217
218 /* Outer loop... */
219 for (it = 1; it <= ret->conv_itmax; it++) {
220
221 /* Store current cost function value... */
222 double chisq_old = chisq;
223
224 /* Compute kernel matrix K_i... */
225 if (it > 1 && it % ret->kernel_recomp == 0)
226 kernel(ctl, tbl, atm_i, obs_i, k_i);
227
228 /* Compute K_i^T * S_eps^{-1} * K_i ... */
229 if (it == 1 || it % ret->kernel_recomp == 0)
230 matrix_product(k_i, sig_eps_inv, 1, cov);
231
232 /* Determine b = K_i^T * S_eps^{-1} * dy - S_a^{-1} * dx ... */
233 for (size_t i = 0; i < m; i++)
234 gsl_vector_set(y_aux, i, gsl_vector_get(dy, i)
235 * POW2(gsl_vector_get(sig_eps_inv, i)));
236 gsl_blas_dgemv(CblasTrans, 1.0, k_i, y_aux, 0.0, b);
237 gsl_blas_dgemv(CblasNoTrans, -1.0, s_a_inv, dx, 1.0, b);
238
239 /* Inner loop... */
240 for (int it2 = 0; it2 < 20; it2++) {
241
242 /* Compute A = (1 + lmpar) * S_a^{-1} + K_i^T * S_eps^{-1} * K_i ... */
243 gsl_matrix_memcpy(a, s_a_inv);
244 gsl_matrix_scale(a, 1 + lmpar);
245 gsl_matrix_add(a, cov);
246
247 /* Solve A * x_step = b by means of Cholesky decomposition... */
248 gsl_linalg_cholesky_decomp(a);
249 gsl_linalg_cholesky_solve(a, b, x_step);
250
251 /* Update atmospheric state... */
252 gsl_vector_add(x_i, x_step);
253 copy_atm(ctl, atm_i, atm_apr, 0);
254 copy_obs(ctl, obs_i, obs_meas, 0);
255 x2atm(ctl, x_i, atm_i);
256
257 /* Check atmospheric state... */
258 for (int ip = 0; ip < atm_i->np; ip++) {
259 atm_i->p[ip] = MIN(MAX(atm_i->p[ip], 5e-7), 5e4);
260 atm_i->t[ip] = MIN(MAX(atm_i->t[ip], 100), 400);
261 for (int ig = 0; ig < ctl->ng; ig++)
262 atm_i->q[ig][ip] = MIN(MAX(atm_i->q[ig][ip], 0), 1);
263 for (int iw = 0; iw < ctl->nw; iw++)
264 atm_i->k[iw][ip] = MAX(atm_i->k[iw][ip], 0);
265 }
266 atm_i->clz = MAX(atm_i->clz, 0);
267 atm_i->cldz = MAX(atm_i->cldz, 0.1);
268 for (int icl = 0; icl < ctl->ncl; icl++)
269 atm_i->clk[icl] = MAX(atm_i->clk[icl], 0);
270 atm_i->sft = MIN(MAX(atm_i->sft, 100), 400);
271 for (int isf = 0; isf < ctl->nsf; isf++)
272 atm_i->sfeps[isf] = MIN(MAX(atm_i->sfeps[isf], 0), 1);
273
274 /* Forward calculation... */
275 formod(ctl, tbl, atm_i, obs_i);
276 obs2y(ctl, obs_i, y_i, NULL, NULL);
277
278 /* Determine dx = x_i - x_a and dy = y - F(x_i) ... */
279 gsl_vector_memcpy(dx, x_i);
280 gsl_vector_sub(dx, x_a);
281 gsl_vector_memcpy(dy, y_m);
282 gsl_vector_sub(dy, y_i);
283
284 /* Compute cost function... */
285 chisq = cost_function(dx, dy, s_a_inv, sig_eps_inv);
286
287 /* Modify Levenberg-Marquardt parameter... */
288 if (chisq > chisq_old) {
289 lmpar *= 10;
290 gsl_vector_sub(x_i, x_step);
291 } else {
292 lmpar /= 10;
293 break;
294 }
295 }
296
297 /* Write info... */
298 LOG(1, "it= %d / chi^2/m= %g", it, chisq);
299
300 /* Write to cost function file... */
301 fprintf(out, "%d %g %d %d\n", it, chisq, (int) m, (int) n);
302
303 /* Get normalized step size in state space... */
304 gsl_blas_ddot(x_step, b, &disq);
305 disq /= (double) n;
306
307 /* Convergence test... */
308 if ((it == 1 || it % ret->kernel_recomp == 0) && disq < ret->conv_dmin)
309 break;
310 }
311
312 /* Close cost function file... */
313 fclose(out);
314
315 /* Store results... */
316 write_atm(ret->dir, "atm_final.tab", ctl, atm_i);
317 write_obs(ret->dir, "obs_final.tab", ctl, obs_i);
318 write_matrix(ret->dir, "matrix_kernel.tab", ctl, k_i,
319 atm_i, obs_i, "y", "x", "r");
320
321 /* ------------------------------------------------------------
322 Analysis of retrieval results...
323 ------------------------------------------------------------ */
324
325 /* Check if error analysis is requested... */
326 if (ret->err_ana) {
327
328 /* Allocate... */
329 gsl_matrix *auxnm = gsl_matrix_alloc(n, m);
330 gsl_matrix *corr = gsl_matrix_alloc(n, n);
331 gsl_matrix *gain = gsl_matrix_alloc(n, m);
332
333 /* Compute inverse retrieval covariance...
334 cov^{-1} = S_a^{-1} + K_i^T * S_eps^{-1} * K_i */
335 matrix_product(k_i, sig_eps_inv, 1, cov);
336 gsl_matrix_add(cov, s_a_inv);
337
338 /* Compute retrieval covariance... */
339 matrix_invert(cov);
340 write_matrix(ret->dir, "matrix_cov_ret.tab", ctl, cov,
341 atm_i, obs_i, "x", "x", "r");
342 write_stddev("total", ret, ctl, atm_i, cov);
343
344 /* Compute correlation matrix... */
345 for (size_t i = 0; i < n; i++)
346 for (size_t j = 0; j < n; j++)
347 gsl_matrix_set(corr, i, j, gsl_matrix_get(cov, i, j)
348 / sqrt(gsl_matrix_get(cov, i, i))
349 / sqrt(gsl_matrix_get(cov, j, j)));
350 write_matrix(ret->dir, "matrix_corr.tab", ctl, corr,
351 atm_i, obs_i, "x", "x", "r");
352
353 /* Compute gain matrix...
354 G = cov * K^T * S_eps^{-1} */
355 for (size_t i = 0; i < n; i++)
356 for (size_t j = 0; j < m; j++)
357 gsl_matrix_set(auxnm, i, j, gsl_matrix_get(k_i, j, i)
358 * POW2(gsl_vector_get(sig_eps_inv, j)));
359 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, cov, auxnm, 0.0, gain);
360 write_matrix(ret->dir, "matrix_gain.tab", ctl, gain,
361 atm_i, obs_i, "x", "y", "c");
362
363 /* Compute retrieval error due to noise... */
364 matrix_product(gain, sig_noise, 2, a);
365 write_stddev("noise", ret, ctl, atm_i, a);
366
367 /* Compute retrieval error due to forward model errors... */
368 matrix_product(gain, sig_formod, 2, a);
369 write_stddev("formod", ret, ctl, atm_i, a);
370
371 /* Compute averaging kernel matrix
372 A = G * K ... */
373 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, gain, k_i, 0.0, a);
374 write_matrix(ret->dir, "matrix_avk.tab", ctl, a,
375 atm_i, obs_i, "x", "x", "r");
376
377 /* Analyze averaging kernel matrix... */
378 analyze_avk(ret, ctl, atm_i, iqa, ipa, a);
379
380 /* Free... */
381 gsl_matrix_free(auxnm);
382 gsl_matrix_free(corr);
383 gsl_matrix_free(gain);
384 }
385
386 /* ------------------------------------------------------------
387 Finalize...
388 ------------------------------------------------------------ */
389
390 gsl_matrix_free(a);
391 gsl_matrix_free(cov);
392 gsl_matrix_free(k_i);
393 gsl_matrix_free(s_a_inv);
394
395 gsl_vector_free(b);
396 gsl_vector_free(dx);
397 gsl_vector_free(dy);
398 gsl_vector_free(sig_eps_inv);
399 gsl_vector_free(sig_formod);
400 gsl_vector_free(sig_noise);
401 gsl_vector_free(x_a);
402 gsl_vector_free(x_i);
403 gsl_vector_free(x_step);
404 gsl_vector_free(y_aux);
405 gsl_vector_free(y_i);
406 gsl_vector_free(y_m);
407}
matrix_product
void matrix_product(gsl_matrix *a, gsl_vector *b, int transpose, gsl_matrix *c)
Definition: jurassic.c:4479
write_atm
void write_atm(const char *dirname, const char *filename, const ctl_t *ctl, const atm_t *atm)
Definition: jurassic.c:5869
set_cov_apr
void set_cov_apr(ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *s_a)
Definition: jurassic.c:5601
x2atm
void x2atm(const ctl_t *ctl, const gsl_vector *x, atm_t *atm)
Definition: jurassic.c:6461
analyze_avk
void analyze_avk(ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *avk)
Definition: jurassic.c:29
set_cov_meas
void set_cov_meas(ret_t *ret, ctl_t *ctl, obs_t *obs, gsl_vector *sig_noise, gsl_vector *sig_formod, gsl_vector *sig_eps_inv)
Definition: jurassic.c:5712
write_stddev
void write_stddev(const char *quantity, ret_t *ret, ctl_t *ctl, atm_t *atm, gsl_matrix *s)
Definition: jurassic.c:6342
copy_obs
void copy_obs(const ctl_t *ctl, obs_t *obs_dest, const obs_t *obs_src, const int init)
Copy or initialize observation geometry and radiance data.
Definition: jurassic.c:3308
write_obs
void write_obs(const char *dirname, const char *filename, const ctl_t *ctl, const obs_t *obs)
Definition: jurassic.c:6202
kernel
void kernel(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute the Jacobian (kernel) matrix by finite differences.
Definition: jurassic.c:4283
matrix_invert
void matrix_invert(gsl_matrix *a)
Definition: jurassic.c:4449
formod
void formod(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs)
Execute the selected forward model.
Definition: jurassic.c:3359
copy_atm
void copy_atm(const ctl_t *ctl, atm_t *atm_dest, const atm_t *atm_src, const int init)
Copy or initialize atmospheric profile data.
Definition: jurassic.c:3258
obs2y
size_t obs2y(const ctl_t *ctl, const obs_t *obs, gsl_vector *y, int *ida, int *ira)
Convert observation radiances into a measurement vector.
Definition: jurassic.c:4524
cost_function
double cost_function(gsl_vector *dx, gsl_vector *dy, gsl_matrix *s_a_inv, gsl_vector *sig_eps_inv)
Definition: jurassic.c:1179
atm2x
size_t atm2x(const ctl_t *ctl, const atm_t *atm, gsl_vector *x, int *iqa, int *ipa)
Convert atmospheric data to state vector elements.
Definition: jurassic.c:110
write_matrix
void write_matrix(const char *dirname, const char *filename, const ctl_t *ctl, const gsl_matrix *matrix, const atm_t *atm, const obs_t *obs, const char *rowspace, const char *colspace, const char *sort)
Definition: jurassic.c:6024
N
#define N
Maximum size of state vector.
Definition: jurassic.h:276
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:266
POW2
#define POW2(x)
Compute the square of a value.
Definition: jurassic.h:669
MIN
#define MIN(a, b)
Determine the minimum of two values.
Definition: jurassic.h:553
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: jurassic.h:864
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: jurassic.h:794
MAX
#define MAX(a, b)
Determine the maximum of two values.
Definition: jurassic.h:535
atm_t
Atmospheric profile data.
Definition: jurassic.h:914
atm_t::sfeps
double sfeps[NSF]
Surface emissivity.
Definition: jurassic.h:956
atm_t::k
double k[NW][NP]
Extinction [km^-1].
Definition: jurassic.h:941
atm_t::t
double t[NP]
Temperature [K].
Definition: jurassic.h:935
atm_t::clz
double clz
Cloud layer height [km].
Definition: jurassic.h:944
atm_t::np
int np
Number of data points.
Definition: jurassic.h:917
atm_t::cldz
double cldz
Cloud layer depth [km].
Definition: jurassic.h:947
atm_t::sft
double sft
Surface temperature [K].
Definition: jurassic.h:953
atm_t::clk
double clk[NCL]
Cloud layer extinction [km^-1].
Definition: jurassic.h:950
atm_t::q
double q[NG][NP]
Volume mixing ratio [ppv].
Definition: jurassic.h:938
atm_t::p
double p[NP]
Pressure [hPa].
Definition: jurassic.h:932
ctl_t::nw
int nw
Number of spectral windows.
Definition: jurassic.h:994
ctl_t::ng
int ng
Number of emitters.
Definition: jurassic.h:970
ctl_t::ncl
int ncl
Number of cloud layer spectral grid points.
Definition: jurassic.h:1000
ctl_t::nsf
int nsf
Number of surface layer spectral grid points.
Definition: jurassic.h:1006
obs_t
Observation geometry and radiance data.
Definition: jurassic.h:1187
ret_t::err_ana
int err_ana
Carry out error analysis (0=no, 1=yes).
Definition: jurassic.h:1254
ret_t::conv_dmin
double conv_dmin
Minimum normalized step size in state space.
Definition: jurassic.h:1251
ret_t::kernel_recomp
int kernel_recomp
Re-computation of kernel matrix (number of iterations).
Definition: jurassic.h:1245
ret_t::conv_itmax
int conv_itmax
Maximum number of iterations.
Definition: jurassic.h:1248
ret_t::dir
char dir[LEN]
Working directory.
Definition: jurassic.h:1242
Here is the call graph for this function:

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 45 of file retrieval.c.

47 {
48
49 static atm_t atm_i, atm_apr;
50 static ctl_t ctl;
51 static obs_t obs_i, obs_meas;
52 static ret_t ret;
53
54 FILE *dirlist;
55
56 /* Check arguments... */
57 if (argc < 3)
58 ERRMSG("Give parameters: <ctl> <dirlist>");
59
60 /* Measure CPU-time... */
61 TIMER("total", 1);
62
63 /* Read control parameters... */
64 read_ctl(argc, argv, &ctl);
65 read_ret(argc, argv, &ctl, &ret);
66
67 /* Initialize look-up tables... */
68 tbl_t *tbl = read_tbl(&ctl);
69
70 /* Open directory list... */
71 if (!(dirlist = fopen(argv[2], "r")))
72 ERRMSG("Cannot open directory list!");
73
74 /* Loop over directories... */
75 while (fscanf(dirlist, "%4999s", ret.dir) != EOF) {
76
77 /* Write info... */
78 LOG(1, "\nRetrieve in directory %s...\n", ret.dir);
79
80 /* Read atmospheric data... */
81 read_atm(ret.dir, "atm_apr.tab", &ctl, &atm_apr);
82
83 /* Read observation data... */
84 read_obs(ret.dir, "obs_meas.tab", &ctl, &obs_meas);
85
86 /* Run retrieval... */
87 optimal_estimation(&ret, &ctl, tbl, &obs_meas, &obs_i, &atm_apr, &atm_i);
88
89 /* Measure CPU-time... */
90 TIMER("total", 2);
91 }
92
93 /* Write info... */
94 LOG(1, "\nRetrieval done...");
95
96 /* Measure CPU-time... */
97 TIMER("total", 3);
98
99 /* Free... */
100 free(tbl);
101
102 return EXIT_SUCCESS;
103}
read_ctl
void read_ctl(int argc, char *argv[], ctl_t *ctl)
Read model control parameters from command-line and configuration input.
Definition: jurassic.c:4889
read_ret
void read_ret(int argc, char *argv[], ctl_t *ctl, ret_t *ret)
Definition: jurassic.c:5202
read_atm
void read_atm(const char *dirname, const char *filename, const ctl_t *ctl, atm_t *atm)
Read atmospheric profile data from an ASCII file.
Definition: jurassic.c:4787
read_obs
void read_obs(const char *dirname, const char *filename, const ctl_t *ctl, obs_t *obs)
Read observation geometry and radiance data from an ASCII file.
Definition: jurassic.c:5043
read_tbl
tbl_t * read_tbl(const ctl_t *ctl)
Read gas emissivity look-up tables for all channels and emitters.
Definition: jurassic.c:5358
TIMER
#define TIMER(name, mode)
Start or stop a named timer.
Definition: jurassic.h:729
optimal_estimation
void optimal_estimation(ret_t *ret, ctl_t *ctl, tbl_t *tbl, obs_t *obs_meas, obs_t *obs_i, atm_t *atm_apr, atm_t *atm_i)
Carry out optimal estimation retrieval.
Definition: retrieval.c:107
ctl_t
Control parameters.
Definition: jurassic.h:967
ret_t
Retrieval control parameters.
Definition: jurassic.h:1239
tbl_t
Emissivity look-up tables.
Definition: jurassic.h:1321
Here is the call graph for this function:
Generated on Wed Nov 12 2025 20:18:22 for JURASSIC by doxygen 1.9.4