JURASSIC
Data Structures | Functions
retrieval.c File Reference

JURASSIC retrieval processor. More...

#include "jurassic.h"

Go to the source code of this file.

Data Structures

struct  ret_t
 Retrieval control parameters. More...
 

Functions

void analyze_avk (ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *avk)
 Compute information content and resolution. More...
 
void analyze_avk_quantity (gsl_matrix *avk, int iq, int *ipa, size_t *n0, size_t *n1, double *cont, double *res)
 Analyze averaging kernels for individual retrieval target. More...
 
double cost_function (gsl_vector *dx, gsl_vector *dy, gsl_matrix *s_a_inv, gsl_vector *sig_eps_inv)
 Compute cost function. More...
 
void matrix_invert (gsl_matrix *a)
 Invert symmetric matrix. More...
 
void matrix_product (gsl_matrix *a, gsl_vector *b, int transpose, gsl_matrix *c)
 Compute matrix product A^TBA or ABA^T for diagonal matrix B. More...
 
void optimal_estimation (ret_t *ret, ctl_t *ctl, obs_t *obs_meas, obs_t *obs_i, atm_t *atm_apr, atm_t *atm_i)
 Carry out optimal estimation retrieval. More...
 
void read_ret (int argc, char *argv[], ctl_t *ctl, ret_t *ret)
 Read retrieval control parameters. More...
 
void set_cov_apr (ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *s_a)
 Set a priori covariance. More...
 
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)
 Set measurement errors. More...
 
void write_stddev (const char *quantity, ret_t *ret, ctl_t *ctl, atm_t *atm, gsl_matrix *s)
 Write retrieval error to file. More...
 
int main (int argc, char *argv[])
 

Detailed Description

JURASSIC retrieval processor.

Definition in file retrieval.c.

Function Documentation

◆ analyze_avk()

void analyze_avk ( ret_t ret,
ctl_t ctl,
atm_t atm,
int *  iqa,
int *  ipa,
gsl_matrix *  avk 
)

Compute information content and resolution.

Definition at line 257 of file retrieval.c.

263 {
264
265 static atm_t atm_cont, atm_res;
266
267 size_t i, n0[NQ], n1[NQ];
268
269 /* Get sizes... */
270 const size_t n = avk->size1;
271
272 /* Find sub-matrices for different quantities... */
273 for (int iq = 0; iq < NQ; iq++) {
274 n0[iq] = N;
275 for (i = 0; i < n; i++) {
276 if (iqa[i] == iq && n0[iq] == N)
277 n0[iq] = i;
278 if (iqa[i] == iq)
279 n1[iq] = i - n0[iq] + 1;
280 }
281 }
282
283 /* Initialize... */
284 copy_atm(ctl, &atm_cont, atm, 1);
285 copy_atm(ctl, &atm_res, atm, 1);
286
287 /* Analyze quantities... */
288 analyze_avk_quantity(avk, IDXP, ipa, n0, n1, atm_cont.p, atm_res.p);
289 analyze_avk_quantity(avk, IDXT, ipa, n0, n1, atm_cont.t, atm_res.t);
290 for (int ig = 0; ig < ctl->ng; ig++)
291 analyze_avk_quantity(avk, IDXQ(ig), ipa, n0, n1,
292 atm_cont.q[ig], atm_res.q[ig]);
293 for (int iw = 0; iw < ctl->nw; iw++)
294 analyze_avk_quantity(avk, IDXK(iw), ipa, n0, n1,
295 atm_cont.k[iw], atm_res.k[iw]);
296 analyze_avk_quantity(avk, IDXCLZ, ipa, n0, n1, &atm_cont.clz, &atm_res.clz);
297 analyze_avk_quantity(avk, IDXCLDZ, ipa, n0, n1, &atm_cont.cldz,
298 &atm_res.cldz);
299 for (int icl = 0; icl < ctl->ncl; icl++)
300 analyze_avk_quantity(avk, IDXCLK(icl), ipa, n0, n1,
301 &atm_cont.clk[icl], &atm_res.clk[icl]);
302 analyze_avk_quantity(avk, IDXSFZ, ipa, n0, n1, &atm_cont.sfz, &atm_res.sfz);
303 analyze_avk_quantity(avk, IDXSFP, ipa, n0, n1, &atm_cont.sfp, &atm_res.sfp);
304 analyze_avk_quantity(avk, IDXSFT, ipa, n0, n1, &atm_cont.sft, &atm_res.sft);
305 for (int isf = 0; isf < ctl->nsf; isf++)
306 analyze_avk_quantity(avk, IDXSFEPS(isf), ipa, n0, n1,
307 &atm_cont.sfeps[isf], &atm_res.sfeps[isf]);
308
309 /* Write results to disk... */
310 write_atm(ret->dir, "atm_cont.tab", ctl, &atm_cont);
311 write_atm(ret->dir, "atm_res.tab", ctl, &atm_res);
312}
void write_atm(const char *dirname, const char *filename, const ctl_t *ctl, const atm_t *atm)
Write atmospheric data.
Definition: jurassic.c:5336
void copy_atm(const ctl_t *ctl, atm_t *atm_dest, const atm_t *atm_src, const int init)
Copy and initialize atmospheric data.
Definition: jurassic.c:2921
#define N
Maximum size of state vector.
Definition: jurassic.h:393
#define IDXCLZ
Index for cloud layer height.
Definition: jurassic.h:463
#define IDXCLDZ
Index for cloud layer depth.
Definition: jurassic.h:466
#define IDXK(iw)
Indices for extinction.
Definition: jurassic.h:460
#define NQ
Maximum number of quantities.
Definition: jurassic.h:398
#define IDXSFT
Index for surface layer temperature.
Definition: jurassic.h:478
#define IDXSFEPS(isf)
Indices for surface layer emissivity.
Definition: jurassic.h:481
#define IDXCLK(icl)
Indices for cloud layer extinction.
Definition: jurassic.h:469
#define IDXP
Index for pressure.
Definition: jurassic.h:451
#define IDXSFP
Index for surface layer pressure.
Definition: jurassic.h:475
#define IDXQ(ig)
Indices for volume mixing ratios.
Definition: jurassic.h:457
#define IDXT
Index for temperature.
Definition: jurassic.h:454
#define IDXSFZ
Index for surface layer height.
Definition: jurassic.h:472
void analyze_avk_quantity(gsl_matrix *avk, int iq, int *ipa, size_t *n0, size_t *n1, double *cont, double *res)
Analyze averaging kernels for individual retrieval target.
Definition: retrieval.c:316
Atmospheric data.
Definition: jurassic.h:488
double sfeps[NSF]
Surface emissivity.
Definition: jurassic.h:536
double k[NW][NP]
Extinction [km^-1].
Definition: jurassic.h:515
double sfz
Surface height [km].
Definition: jurassic.h:527
double t[NP]
Temperature [K].
Definition: jurassic.h:509
double sfp
Surface pressure [hPa].
Definition: jurassic.h:530
double clz
Cloud layer height [km].
Definition: jurassic.h:518
double cldz
Cloud layer depth [km].
Definition: jurassic.h:521
double sft
Surface temperature [K].
Definition: jurassic.h:533
double clk[NCL]
Cloud layer extinction [km^-1].
Definition: jurassic.h:524
double q[NG][NP]
Volume mixing ratio [ppv].
Definition: jurassic.h:512
double p[NP]
Pressure [hPa].
Definition: jurassic.h:506
int nw
Number of spectral windows.
Definition: jurassic.h:556
int ng
Number of emitters.
Definition: jurassic.h:544
int ncl
Number of cloud layer spectral grid points.
Definition: jurassic.h:562
int nsf
Number of surface layer spectral grid points.
Definition: jurassic.h:568
char dir[LEN]
Working directory.
Definition: retrieval.c:35
Here is the call graph for this function:

◆ analyze_avk_quantity()

void analyze_avk_quantity ( gsl_matrix *  avk,
int  iq,
int *  ipa,
size_t *  n0,
size_t *  n1,
double *  cont,
double *  res 
)

Analyze averaging kernels for individual retrieval target.

Definition at line 316 of file retrieval.c.

323 {
324
325 /* Loop over state vector elements... */
326 if (n0[iq] < N)
327 for (size_t i = 0; i < n1[iq]; i++) {
328
329 /* Get area of averaging kernel... */
330 for (size_t j = 0; j < n1[iq]; j++)
331 cont[ipa[n0[iq] + i]] += gsl_matrix_get(avk, n0[iq] + i, n0[iq] + j);
332
333 /* Get information density... */
334 res[ipa[n0[iq] + i]] = 1 / gsl_matrix_get(avk, n0[iq] + i, n0[iq] + i);
335 }
336}

◆ cost_function()

double cost_function ( gsl_vector *  dx,
gsl_vector *  dy,
gsl_matrix *  s_a_inv,
gsl_vector *  sig_eps_inv 
)

Compute cost function.

Definition at line 340 of file retrieval.c.

344 {
345
346 double chisq_a, chisq_m = 0;
347
348 /* Get sizes... */
349 const size_t m = dy->size;
350 const size_t n = dx->size;
351
352 /* Allocate... */
353 gsl_vector *x_aux = gsl_vector_alloc(n);
354 gsl_vector *y_aux = gsl_vector_alloc(m);
355
356 /* Determine normalized cost function...
357 (chi^2 = 1/m * [dy^T * S_eps^{-1} * dy + dx^T * S_a^{-1} * dx]) */
358 for (size_t i = 0; i < m; i++)
359 chisq_m += POW2(gsl_vector_get(dy, i) * gsl_vector_get(sig_eps_inv, i));
360 gsl_blas_dgemv(CblasNoTrans, 1.0, s_a_inv, dx, 0.0, x_aux);
361 gsl_blas_ddot(dx, x_aux, &chisq_a);
362
363 /* Free... */
364 gsl_vector_free(x_aux);
365 gsl_vector_free(y_aux);
366
367 /* Return cost function value... */
368 return (chisq_m + chisq_a) / (double) m;
369}
#define POW2(x)
Compute x^2.
Definition: jurassic.h:187

◆ matrix_invert()

void matrix_invert ( gsl_matrix *  a)

Invert symmetric matrix.

Definition at line 373 of file retrieval.c.

374 {
375
376 size_t diag = 1;
377
378 /* Get size... */
379 const size_t n = a->size1;
380
381 /* Check if matrix is diagonal... */
382 for (size_t i = 0; i < n && diag; i++)
383 for (size_t j = i + 1; j < n; j++)
384 if (gsl_matrix_get(a, i, j) != 0) {
385 diag = 0;
386 break;
387 }
388
389 /* Quick inversion of diagonal matrix... */
390 if (diag)
391 for (size_t i = 0; i < n; i++)
392 gsl_matrix_set(a, i, i, 1 / gsl_matrix_get(a, i, i));
393
394 /* Matrix inversion by means of Cholesky decomposition... */
395 else {
396 gsl_linalg_cholesky_decomp(a);
397 gsl_linalg_cholesky_invert(a);
398 }
399}

◆ matrix_product()

void matrix_product ( gsl_matrix *  a,
gsl_vector *  b,
int  transpose,
gsl_matrix *  c 
)

Compute matrix product A^TBA or ABA^T for diagonal matrix B.

Definition at line 403 of file retrieval.c.

407 {
408
409 /* Set sizes... */
410 const size_t m = a->size1;
411 const size_t n = a->size2;
412
413 /* Allocate... */
414 gsl_matrix *aux = gsl_matrix_alloc(m, n);
415
416 /* Compute A^T B A... */
417 if (transpose == 1) {
418
419 /* Compute B^1/2 A... */
420 for (size_t i = 0; i < m; i++)
421 for (size_t j = 0; j < n; j++)
422 gsl_matrix_set(aux, i, j,
423 gsl_vector_get(b, i) * gsl_matrix_get(a, i, j));
424
425 /* Compute A^T B A = (B^1/2 A)^T (B^1/2 A)... */
426 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, aux, aux, 0.0, c);
427 }
428
429 /* Compute A B A^T... */
430 else if (transpose == 2) {
431
432 /* Compute A B^1/2... */
433 for (size_t i = 0; i < m; i++)
434 for (size_t j = 0; j < n; j++)
435 gsl_matrix_set(aux, i, j,
436 gsl_matrix_get(a, i, j) * gsl_vector_get(b, j));
437
438 /* Compute A B A^T = (A B^1/2) (A B^1/2)^T... */
439 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, aux, aux, 0.0, c);
440 }
441
442 /* Free... */
443 gsl_matrix_free(aux);
444}

◆ optimal_estimation()

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

Carry out optimal estimation retrieval.

Definition at line 448 of file retrieval.c.

454 {
455
456 static int ipa[N], iqa[N];
457
458 FILE *out;
459
460 char filename[2 * LEN];
461
462 double chisq, disq = 0, lmpar = 0.001;
463
464 int it = 0;
465
466 /* ------------------------------------------------------------
467 Initialize...
468 ------------------------------------------------------------ */
469
470 /* Get sizes... */
471 const size_t m = obs2y(ctl, obs_meas, NULL, NULL, NULL);
472 const size_t n = atm2x(ctl, atm_apr, NULL, iqa, ipa);
473 if (m == 0 || n == 0)
474 ERRMSG("Check problem definition!");
475
476 /* Write info... */
477 LOG(1, "Problem size: m= %d / n= %d "
478 "(alloc= %.4g MB / stat= %.4g MB)",
479 (int) m, (int) n,
480 (double) (3 * m * n + 4 * n * n + 8 * m +
481 8 * n) * sizeof(double) / 1024. / 1024.,
482 (double) (5 * sizeof(atm_t) + 3 * sizeof(obs_t)
483 + 2 * N * sizeof(int)) / 1024. / 1024.);
484
485 /* Allocate... */
486 gsl_matrix *a = gsl_matrix_alloc(n, n);
487 gsl_matrix *cov = gsl_matrix_alloc(n, n);
488 gsl_matrix *k_i = gsl_matrix_alloc(m, n);
489 gsl_matrix *s_a_inv = gsl_matrix_alloc(n, n);
490
491 gsl_vector *b = gsl_vector_alloc(n);
492 gsl_vector *dx = gsl_vector_alloc(n);
493 gsl_vector *dy = gsl_vector_alloc(m);
494 gsl_vector *sig_eps_inv = gsl_vector_alloc(m);
495 gsl_vector *sig_formod = gsl_vector_alloc(m);
496 gsl_vector *sig_noise = gsl_vector_alloc(m);
497 gsl_vector *x_a = gsl_vector_alloc(n);
498 gsl_vector *x_i = gsl_vector_alloc(n);
499 gsl_vector *x_step = gsl_vector_alloc(n);
500 gsl_vector *y_aux = gsl_vector_alloc(m);
501 gsl_vector *y_i = gsl_vector_alloc(m);
502 gsl_vector *y_m = gsl_vector_alloc(m);
503
504 /* Set initial state... */
505 copy_atm(ctl, atm_i, atm_apr, 0);
506 copy_obs(ctl, obs_i, obs_meas, 0);
507 formod(ctl, atm_i, obs_i);
508
509 /* Set state vectors and observation vectors... */
510 atm2x(ctl, atm_apr, x_a, NULL, NULL);
511 atm2x(ctl, atm_i, x_i, NULL, NULL);
512 obs2y(ctl, obs_meas, y_m, NULL, NULL);
513 obs2y(ctl, obs_i, y_i, NULL, NULL);
514
515 /* Set inverse a priori covariance S_a^-1... */
516 set_cov_apr(ret, ctl, atm_apr, iqa, ipa, s_a_inv);
517 write_matrix(ret->dir, "matrix_cov_apr.tab", ctl, s_a_inv,
518 atm_i, obs_i, "x", "x", "r");
519 matrix_invert(s_a_inv);
520
521 /* Get measurement errors... */
522 set_cov_meas(ret, ctl, obs_meas, sig_noise, sig_formod, sig_eps_inv);
523
524 /* Create cost function file... */
525 sprintf(filename, "%s/costs.tab", ret->dir);
526 if (!(out = fopen(filename, "w")))
527 ERRMSG("Cannot create cost function file!");
528
529 /* Write header... */
530 fprintf(out,
531 "# $1 = iteration number\n"
532 "# $2 = normalized cost function\n"
533 "# $3 = number of measurements\n"
534 "# $4 = number of state vector elements\n\n");
535
536 /* Determine dx = x_i - x_a and dy = y - F(x_i) ... */
537 gsl_vector_memcpy(dx, x_i);
538 gsl_vector_sub(dx, x_a);
539 gsl_vector_memcpy(dy, y_m);
540 gsl_vector_sub(dy, y_i);
541
542 /* Compute cost function... */
543 chisq = cost_function(dx, dy, s_a_inv, sig_eps_inv);
544
545 /* Write info... */
546 LOG(1, "it= %d / chi^2/m= %g", it, chisq);
547
548 /* Write to cost function file... */
549 fprintf(out, "%d %g %d %d\n", it, chisq, (int) m, (int) n);
550
551 /* Compute initial kernel... */
552 kernel(ctl, atm_i, obs_i, k_i);
553
554 /* ------------------------------------------------------------
555 Levenberg-Marquardt minimization...
556 ------------------------------------------------------------ */
557
558 /* Outer loop... */
559 for (it = 1; it <= ret->conv_itmax; it++) {
560
561 /* Store current cost function value... */
562 double chisq_old = chisq;
563
564 /* Compute kernel matrix K_i... */
565 if (it > 1 && it % ret->kernel_recomp == 0)
566 kernel(ctl, atm_i, obs_i, k_i);
567
568 /* Compute K_i^T * S_eps^{-1} * K_i ... */
569 if (it == 1 || it % ret->kernel_recomp == 0)
570 matrix_product(k_i, sig_eps_inv, 1, cov);
571
572 /* Determine b = K_i^T * S_eps^{-1} * dy - S_a^{-1} * dx ... */
573 for (size_t i = 0; i < m; i++)
574 gsl_vector_set(y_aux, i, gsl_vector_get(dy, i)
575 * POW2(gsl_vector_get(sig_eps_inv, i)));
576 gsl_blas_dgemv(CblasTrans, 1.0, k_i, y_aux, 0.0, b);
577 gsl_blas_dgemv(CblasNoTrans, -1.0, s_a_inv, dx, 1.0, b);
578
579 /* Inner loop... */
580 for (int it2 = 0; it2 < 20; it2++) {
581
582 /* Compute A = (1 + lmpar) * S_a^{-1} + K_i^T * S_eps^{-1} * K_i ... */
583 gsl_matrix_memcpy(a, s_a_inv);
584 gsl_matrix_scale(a, 1 + lmpar);
585 gsl_matrix_add(a, cov);
586
587 /* Solve A * x_step = b by means of Cholesky decomposition... */
588 gsl_linalg_cholesky_decomp(a);
589 gsl_linalg_cholesky_solve(a, b, x_step);
590
591 /* Update atmospheric state... */
592 gsl_vector_add(x_i, x_step);
593 copy_atm(ctl, atm_i, atm_apr, 0);
594 copy_obs(ctl, obs_i, obs_meas, 0);
595 x2atm(ctl, x_i, atm_i);
596
597 /* Check atmospheric state... */
598 for (int ip = 0; ip < atm_i->np; ip++) {
599 atm_i->p[ip] = MIN(MAX(atm_i->p[ip], 5e-7), 5e4);
600 atm_i->t[ip] = MIN(MAX(atm_i->t[ip], 100), 400);
601 for (int ig = 0; ig < ctl->ng; ig++)
602 atm_i->q[ig][ip] = MIN(MAX(atm_i->q[ig][ip], 0), 1);
603 for (int iw = 0; iw < ctl->nw; iw++)
604 atm_i->k[iw][ip] = MAX(atm_i->k[iw][ip], 0);
605 }
606 atm_i->clz = MAX(atm_i->clz, 0);
607 atm_i->cldz = MAX(atm_i->cldz, 0.1);
608 for (int icl = 0; icl < ctl->ncl; icl++)
609 atm_i->clk[icl] = MAX(atm_i->clk[icl], 0);
610 atm_i->sfz = MAX(atm_i->sfz, 0);
611 atm_i->sfp = MAX(atm_i->sfp, 0);
612 atm_i->sft = MIN(MAX(atm_i->sft, 100), 400);
613 for (int isf = 0; isf < ctl->nsf; isf++)
614 atm_i->sfeps[isf] = MIN(MAX(atm_i->sfeps[isf], 0), 1);
615
616 /* Forward calculation... */
617 formod(ctl, atm_i, obs_i);
618 obs2y(ctl, obs_i, y_i, NULL, NULL);
619
620 /* Determine dx = x_i - x_a and dy = y - F(x_i) ... */
621 gsl_vector_memcpy(dx, x_i);
622 gsl_vector_sub(dx, x_a);
623 gsl_vector_memcpy(dy, y_m);
624 gsl_vector_sub(dy, y_i);
625
626 /* Compute cost function... */
627 chisq = cost_function(dx, dy, s_a_inv, sig_eps_inv);
628
629 /* Modify Levenberg-Marquardt parameter... */
630 if (chisq > chisq_old) {
631 lmpar *= 10;
632 gsl_vector_sub(x_i, x_step);
633 } else {
634 lmpar /= 10;
635 break;
636 }
637 }
638
639 /* Write info... */
640 LOG(1, "it= %d / chi^2/m= %g", it, chisq);
641
642 /* Write to cost function file... */
643 fprintf(out, "%d %g %d %d\n", it, chisq, (int) m, (int) n);
644
645 /* Get normalized step size in state space... */
646 gsl_blas_ddot(x_step, b, &disq);
647 disq /= (double) n;
648
649 /* Convergence test... */
650 if ((it == 1 || it % ret->kernel_recomp == 0) && disq < ret->conv_dmin)
651 break;
652 }
653
654 /* Close cost function file... */
655 fclose(out);
656
657 /* Store results... */
658 write_atm(ret->dir, "atm_final.tab", ctl, atm_i);
659 write_obs(ret->dir, "obs_final.tab", ctl, obs_i);
660 write_matrix(ret->dir, "matrix_kernel.tab", ctl, k_i,
661 atm_i, obs_i, "y", "x", "r");
662
663 /* ------------------------------------------------------------
664 Analysis of retrieval results...
665 ------------------------------------------------------------ */
666
667 /* Check if error analysis is requested... */
668 if (ret->err_ana) {
669
670 /* Allocate... */
671 gsl_matrix *auxnm = gsl_matrix_alloc(n, m);
672 gsl_matrix *corr = gsl_matrix_alloc(n, n);
673 gsl_matrix *gain = gsl_matrix_alloc(n, m);
674
675 /* Compute inverse retrieval covariance...
676 cov^{-1} = S_a^{-1} + K_i^T * S_eps^{-1} * K_i */
677 matrix_product(k_i, sig_eps_inv, 1, cov);
678 gsl_matrix_add(cov, s_a_inv);
679
680 /* Compute retrieval covariance... */
681 matrix_invert(cov);
682 write_matrix(ret->dir, "matrix_cov_ret.tab", ctl, cov,
683 atm_i, obs_i, "x", "x", "r");
684 write_stddev("total", ret, ctl, atm_i, cov);
685
686 /* Compute correlation matrix... */
687 for (size_t i = 0; i < n; i++)
688 for (size_t j = 0; j < n; j++)
689 gsl_matrix_set(corr, i, j, gsl_matrix_get(cov, i, j)
690 / sqrt(gsl_matrix_get(cov, i, i))
691 / sqrt(gsl_matrix_get(cov, j, j)));
692 write_matrix(ret->dir, "matrix_corr.tab", ctl, corr,
693 atm_i, obs_i, "x", "x", "r");
694
695 /* Compute gain matrix...
696 G = cov * K^T * S_eps^{-1} */
697 for (size_t i = 0; i < n; i++)
698 for (size_t j = 0; j < m; j++)
699 gsl_matrix_set(auxnm, i, j, gsl_matrix_get(k_i, j, i)
700 * POW2(gsl_vector_get(sig_eps_inv, j)));
701 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, cov, auxnm, 0.0, gain);
702 write_matrix(ret->dir, "matrix_gain.tab", ctl, gain,
703 atm_i, obs_i, "x", "y", "c");
704
705 /* Compute retrieval error due to noise... */
706 matrix_product(gain, sig_noise, 2, a);
707 write_stddev("noise", ret, ctl, atm_i, a);
708
709 /* Compute retrieval error due to forward model errors... */
710 matrix_product(gain, sig_formod, 2, a);
711 write_stddev("formod", ret, ctl, atm_i, a);
712
713 /* Compute averaging kernel matrix
714 A = G * K ... */
715 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, gain, k_i, 0.0, a);
716 write_matrix(ret->dir, "matrix_avk.tab", ctl, a,
717 atm_i, obs_i, "x", "x", "r");
718
719 /* Analyze averaging kernel matrix... */
720 analyze_avk(ret, ctl, atm_i, iqa, ipa, a);
721
722 /* Free... */
723 gsl_matrix_free(auxnm);
724 gsl_matrix_free(corr);
725 gsl_matrix_free(gain);
726 }
727
728 /* ------------------------------------------------------------
729 Finalize...
730 ------------------------------------------------------------ */
731
732 gsl_matrix_free(a);
733 gsl_matrix_free(cov);
734 gsl_matrix_free(k_i);
735 gsl_matrix_free(s_a_inv);
736
737 gsl_vector_free(b);
738 gsl_vector_free(dx);
739 gsl_vector_free(dy);
740 gsl_vector_free(sig_eps_inv);
741 gsl_vector_free(sig_formod);
742 gsl_vector_free(sig_noise);
743 gsl_vector_free(x_a);
744 gsl_vector_free(x_i);
745 gsl_vector_free(x_step);
746 gsl_vector_free(y_aux);
747 gsl_vector_free(y_i);
748 gsl_vector_free(y_m);
749}
void x2atm(const ctl_t *ctl, const gsl_vector *x, atm_t *atm)
Decompose parameter vector or state vector.
Definition: jurassic.c:5896
void formod(const ctl_t *ctl, atm_t *atm, obs_t *obs)
Determine ray paths and compute radiative transfer.
Definition: jurassic.c:3026
void copy_obs(const ctl_t *ctl, obs_t *obs_dest, const obs_t *obs_src, const int init)
Copy and initialize observation data.
Definition: jurassic.c:2975
void write_obs(const char *dirname, const char *filename, const ctl_t *ctl, const obs_t *obs)
Write observation data.
Definition: jurassic.c:5670
size_t obs2y(const ctl_t *ctl, const obs_t *obs, gsl_vector *y, int *ida, int *ira)
Compose measurement vector.
Definition: jurassic.c:4174
size_t atm2x(const ctl_t *ctl, const atm_t *atm, gsl_vector *x, int *iqa, int *ipa)
Compose state vector or parameter vector.
Definition: jurassic.c:29
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)
Write matrix.
Definition: jurassic.c:5492
void kernel(ctl_t *ctl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute Jacobians.
Definition: jurassic.c:4004
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:383
#define MIN(a, b)
Macro to determine the minimum of two values.
Definition: jurassic.h:160
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
#define LOG(level,...)
Print log message.
Definition: jurassic.h:221
#define MAX(a, b)
Macro to determine the maximum of two values.
Definition: jurassic.h:156
void matrix_product(gsl_matrix *a, gsl_vector *b, int transpose, gsl_matrix *c)
Compute matrix product A^TBA or ABA^T for diagonal matrix B.
Definition: retrieval.c:403
void set_cov_apr(ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *s_a)
Set a priori covariance.
Definition: retrieval.c:808
void analyze_avk(ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *avk)
Compute information content and resolution.
Definition: retrieval.c:257
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)
Set measurement errors.
Definition: retrieval.c:924
void write_stddev(const char *quantity, ret_t *ret, ctl_t *ctl, atm_t *atm, gsl_matrix *s)
Write retrieval error to file.
Definition: retrieval.c:968
void matrix_invert(gsl_matrix *a)
Invert symmetric matrix.
Definition: retrieval.c:373
double cost_function(gsl_vector *dx, gsl_vector *dy, gsl_matrix *s_a_inv, gsl_vector *sig_eps_inv)
Compute cost function.
Definition: retrieval.c:340
int np
Number of data points.
Definition: jurassic.h:491
Observation geometry and radiance data.
Definition: jurassic.h:734
int err_ana
Carry out error analysis (0=no, 1=yes).
Definition: retrieval.c:47
double conv_dmin
Minimum normalized step size in state space.
Definition: retrieval.c:44
int kernel_recomp
Re-computation of kernel matrix (number of iterations).
Definition: retrieval.c:38
int conv_itmax
Maximum number of iterations.
Definition: retrieval.c:41
Here is the call graph for this function:

◆ read_ret()

void read_ret ( int  argc,
char *  argv[],
ctl_t ctl,
ret_t ret 
)

Read retrieval control parameters.

Definition at line 753 of file retrieval.c.

757 {
758
759 /* Iteration control... */
760 ret->kernel_recomp =
761 (int) scan_ctl(argc, argv, "KERNEL_RECOMP", -1, "3", NULL);
762 ret->conv_itmax = (int) scan_ctl(argc, argv, "CONV_ITMAX", -1, "30", NULL);
763 ret->conv_dmin = scan_ctl(argc, argv, "CONV_DMIN", -1, "0.1", NULL);
764
765 /* Error analysis... */
766 ret->err_ana = (int) scan_ctl(argc, argv, "ERR_ANA", -1, "1", NULL);
767
768 for (int id = 0; id < ctl->nd; id++)
769 ret->err_formod[id] = scan_ctl(argc, argv, "ERR_FORMOD", id, "0", NULL);
770
771 for (int id = 0; id < ctl->nd; id++)
772 ret->err_noise[id] = scan_ctl(argc, argv, "ERR_NOISE", id, "0", NULL);
773
774 ret->err_press = scan_ctl(argc, argv, "ERR_PRESS", -1, "0", NULL);
775 ret->err_press_cz = scan_ctl(argc, argv, "ERR_PRESS_CZ", -1, "-999", NULL);
776 ret->err_press_ch = scan_ctl(argc, argv, "ERR_PRESS_CH", -1, "-999", NULL);
777
778 ret->err_temp = scan_ctl(argc, argv, "ERR_TEMP", -1, "0", NULL);
779 ret->err_temp_cz = scan_ctl(argc, argv, "ERR_TEMP_CZ", -1, "-999", NULL);
780 ret->err_temp_ch = scan_ctl(argc, argv, "ERR_TEMP_CH", -1, "-999", NULL);
781
782 for (int ig = 0; ig < ctl->ng; ig++) {
783 ret->err_q[ig] = scan_ctl(argc, argv, "ERR_Q", ig, "0", NULL);
784 ret->err_q_cz[ig] = scan_ctl(argc, argv, "ERR_Q_CZ", ig, "-999", NULL);
785 ret->err_q_ch[ig] = scan_ctl(argc, argv, "ERR_Q_CH", ig, "-999", NULL);
786 }
787
788 for (int iw = 0; iw < ctl->nw; iw++) {
789 ret->err_k[iw] = scan_ctl(argc, argv, "ERR_K", iw, "0", NULL);
790 ret->err_k_cz[iw] = scan_ctl(argc, argv, "ERR_K_CZ", iw, "-999", NULL);
791 ret->err_k_ch[iw] = scan_ctl(argc, argv, "ERR_K_CH", iw, "-999", NULL);
792 }
793
794 ret->err_clz = scan_ctl(argc, argv, "ERR_CLZ", -1, "0", NULL);
795 ret->err_cldz = scan_ctl(argc, argv, "ERR_CLDZ", -1, "0", NULL);
796 for (int icl = 0; icl < ctl->ncl; icl++)
797 ret->err_clk[icl] = scan_ctl(argc, argv, "ERR_CLK", icl, "0", NULL);
798
799 ret->err_sfz = scan_ctl(argc, argv, "ERR_SFZ", -1, "0", NULL);
800 ret->err_sfp = scan_ctl(argc, argv, "ERR_SFP", -1, "0", NULL);
801 ret->err_sft = scan_ctl(argc, argv, "ERR_SFT", -1, "0", NULL);
802 for (int isf = 0; isf < ctl->nsf; isf++)
803 ret->err_sfeps[isf] = scan_ctl(argc, argv, "ERR_SFEPS", isf, "0", NULL);
804}
double scan_ctl(int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
Definition: jurassic.c:5114
int nd
Number of radiance channels.
Definition: jurassic.h:550
double err_sfz
Surface height error [km].
Definition: retrieval.c:101
double err_press_cz
Vertical correlation length for pressure error [km].
Definition: retrieval.c:59
double err_sfp
Surface pressure error [hPa].
Definition: retrieval.c:104
double err_press
Pressure error [%].
Definition: retrieval.c:56
double err_k_cz[NW]
Vertical correlation length for extinction error [km].
Definition: retrieval.c:86
double err_k_ch[NW]
Horizontal correlation length for extinction error [km].
Definition: retrieval.c:89
double err_temp_cz
Vertical correlation length for temperature error [km].
Definition: retrieval.c:68
double err_formod[ND]
Forward model error [%].
Definition: retrieval.c:50
double err_temp
Temperature error [K].
Definition: retrieval.c:65
double err_q_cz[NG]
Vertical correlation length for volume mixing ratio error [km].
Definition: retrieval.c:77
double err_noise[ND]
Noise error [W/(m^2 sr cm^-1)].
Definition: retrieval.c:53
double err_clz
Cloud height error [km].
Definition: retrieval.c:92
double err_sft
Surface temperature error [K].
Definition: retrieval.c:107
double err_clk[NCL]
Cloud extinction error [km^-1].
Definition: retrieval.c:98
double err_temp_ch
Horizontal correlation length for temperature error [km].
Definition: retrieval.c:71
double err_cldz
Cloud depth error [km].
Definition: retrieval.c:95
double err_press_ch
Horizontal correlation length for pressure error [km].
Definition: retrieval.c:62
double err_q_ch[NG]
Horizontal correlation length for volume mixing ratio error [km].
Definition: retrieval.c:80
double err_q[NG]
Volume mixing ratio error [%].
Definition: retrieval.c:74
double err_sfeps[NSF]
Surface emissivity error.
Definition: retrieval.c:110
double err_k[NW]
Extinction error [km^-1].
Definition: retrieval.c:83
Here is the call graph for this function:

◆ set_cov_apr()

void set_cov_apr ( ret_t ret,
ctl_t ctl,
atm_t atm,
int *  iqa,
int *  ipa,
gsl_matrix *  s_a 
)

Set a priori covariance.

Definition at line 808 of file retrieval.c.

814 {
815
816 double x0[3], x1[3];
817
818 /* Get sizes... */
819 const size_t n = s_a->size1;
820
821 /* Allocate... */
822 gsl_vector *x_a = gsl_vector_alloc(n);
823
824 /* Get sigma vector... */
825 atm2x(ctl, atm, x_a, NULL, NULL);
826 for (size_t i = 0; i < n; i++) {
827 if (iqa[i] == IDXP)
828 gsl_vector_set(x_a, i, ret->err_press / 100 * gsl_vector_get(x_a, i));
829 if (iqa[i] == IDXT)
830 gsl_vector_set(x_a, i, ret->err_temp);
831 for (int ig = 0; ig < ctl->ng; ig++)
832 if (iqa[i] == IDXQ(ig))
833 gsl_vector_set(x_a, i, ret->err_q[ig] / 100 * gsl_vector_get(x_a, i));
834 for (int iw = 0; iw < ctl->nw; iw++)
835 if (iqa[i] == IDXK(iw))
836 gsl_vector_set(x_a, i, ret->err_k[iw]);
837 if (iqa[i] == IDXCLZ)
838 gsl_vector_set(x_a, i, ret->err_clz);
839 if (iqa[i] == IDXCLDZ)
840 gsl_vector_set(x_a, i, ret->err_cldz);
841 for (int icl = 0; icl < ctl->ncl; icl++)
842 if (iqa[i] == IDXCLK(icl))
843 gsl_vector_set(x_a, i, ret->err_clk[icl]);
844 if (iqa[i] == IDXSFZ)
845 gsl_vector_set(x_a, i, ret->err_sfz);
846 if (iqa[i] == IDXSFP)
847 gsl_vector_set(x_a, i, ret->err_sfp);
848 if (iqa[i] == IDXSFT)
849 gsl_vector_set(x_a, i, ret->err_sft);
850 for (int isf = 0; isf < ctl->nsf; isf++)
851 if (iqa[i] == IDXSFEPS(isf))
852 gsl_vector_set(x_a, i, ret->err_sfeps[isf]);
853 }
854
855 /* Check standard deviations... */
856 for (size_t i = 0; i < n; i++)
857 if (POW2(gsl_vector_get(x_a, i)) <= 0)
858 ERRMSG("Check a priori data (zero standard deviation)!");
859
860 /* Initialize diagonal covariance... */
861 gsl_matrix_set_zero(s_a);
862 for (size_t i = 0; i < n; i++)
863 gsl_matrix_set(s_a, i, i, POW2(gsl_vector_get(x_a, i)));
864
865 /* Loop over matrix elements... */
866 for (size_t i = 0; i < n; i++)
867 for (size_t j = 0; j < n; j++)
868 if (i != j && iqa[i] == iqa[j]) {
869
870 /* Initialize... */
871 double cz = 0;
872 double ch = 0;
873
874 /* Set correlation lengths for pressure... */
875 if (iqa[i] == IDXP) {
876 cz = ret->err_press_cz;
877 ch = ret->err_press_ch;
878 }
879
880 /* Set correlation lengths for temperature... */
881 if (iqa[i] == IDXT) {
882 cz = ret->err_temp_cz;
883 ch = ret->err_temp_ch;
884 }
885
886 /* Set correlation lengths for volume mixing ratios... */
887 for (int ig = 0; ig < ctl->ng; ig++)
888 if (iqa[i] == IDXQ(ig)) {
889 cz = ret->err_q_cz[ig];
890 ch = ret->err_q_ch[ig];
891 }
892
893 /* Set correlation lengths for extinction... */
894 for (int iw = 0; iw < ctl->nw; iw++)
895 if (iqa[i] == IDXK(iw)) {
896 cz = ret->err_k_cz[iw];
897 ch = ret->err_k_ch[iw];
898 }
899
900 /* Compute correlations... */
901 if (cz > 0 && ch > 0) {
902
903 /* Get Cartesian coordinates... */
904 geo2cart(0, atm->lon[ipa[i]], atm->lat[ipa[i]], x0);
905 geo2cart(0, atm->lon[ipa[j]], atm->lat[ipa[j]], x1);
906
907 /* Compute correlations... */
908 double rho =
909 exp(-DIST(x0, x1) / ch -
910 fabs(atm->z[ipa[i]] - atm->z[ipa[j]]) / cz);
911
912 /* Set covariance... */
913 gsl_matrix_set(s_a, i, j, gsl_vector_get(x_a, i)
914 * gsl_vector_get(x_a, j) * rho);
915 }
916 }
917
918 /* Free... */
919 gsl_vector_free(x_a);
920}
void geo2cart(const double z, const double lon, const double lat, double *x)
Convert geolocation to Cartesian coordinates.
Definition: jurassic.c:3500
#define DIST(a, b)
Compute Cartesian distance between two vectors.
Definition: jurassic.h:134
double lat[NP]
Latitude [deg].
Definition: jurassic.h:503
double lon[NP]
Longitude [deg].
Definition: jurassic.h:500
double z[NP]
Altitude [km].
Definition: jurassic.h:497
Here is the call graph for this function:

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

Set measurement errors.

Definition at line 924 of file retrieval.c.

930 {
931
932 static obs_t obs_err;
933
934 /* Get size... */
935 const size_t m = sig_eps_inv->size;
936
937 /* Noise error (always considered in retrieval fit)... */
938 copy_obs(ctl, &obs_err, obs, 1);
939 for (int ir = 0; ir < obs_err.nr; ir++)
940 for (int id = 0; id < ctl->nd; id++)
941 obs_err.rad[id][ir]
942 = (isfinite(obs->rad[id][ir]) ? ret->err_noise[id] : NAN);
943 obs2y(ctl, &obs_err, sig_noise, NULL, NULL);
944
945 /* Forward model error (always considered in retrieval fit)... */
946 copy_obs(ctl, &obs_err, obs, 1);
947 for (int ir = 0; ir < obs_err.nr; ir++)
948 for (int id = 0; id < ctl->nd; id++)
949 obs_err.rad[id][ir]
950 = fabs(ret->err_formod[id] / 100 * obs->rad[id][ir]);
951 obs2y(ctl, &obs_err, sig_formod, NULL, NULL);
952
953 /* Total error... */
954 for (size_t i = 0; i < m; i++)
955 gsl_vector_set(sig_eps_inv, i, 1 / sqrt(POW2(gsl_vector_get(sig_noise, i))
956 +
957 POW2(gsl_vector_get
958 (sig_formod, i))));
959
960 /* Check standard deviations... */
961 for (size_t i = 0; i < m; i++)
962 if (gsl_vector_get(sig_eps_inv, i) <= 0)
963 ERRMSG("Check measurement errors (zero standard deviation)!");
964}
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
Definition: jurassic.h:773
int nr
Number of ray paths.
Definition: jurassic.h:737
Here is the call graph for this function:

◆ write_stddev()

void write_stddev ( const char *  quantity,
ret_t ret,
ctl_t ctl,
atm_t atm,
gsl_matrix *  s 
)

Write retrieval error to file.

Definition at line 968 of file retrieval.c.

973 {
974
975 static atm_t atm_aux;
976
977 char filename[LEN];
978
979 /* Get sizes... */
980 const size_t n = s->size1;
981
982 /* Allocate... */
983 gsl_vector *x_aux = gsl_vector_alloc(n);
984
985 /* Compute standard deviation... */
986 for (size_t i = 0; i < n; i++)
987 gsl_vector_set(x_aux, i, sqrt(gsl_matrix_get(s, i, i)));
988
989 /* Write to disk... */
990 copy_atm(ctl, &atm_aux, atm, 1);
991 x2atm(ctl, x_aux, &atm_aux);
992 sprintf(filename, "atm_err_%s.tab", quantity);
993 write_atm(ret->dir, filename, ctl, &atm_aux);
994
995 /* Free... */
996 gsl_vector_free(x_aux);
997}
Here is the call graph for this function:

◆ main()

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

Definition at line 201 of file retrieval.c.

203 {
204
205 static atm_t atm_i, atm_apr;
206 static ctl_t ctl;
207 static obs_t obs_i, obs_meas;
208 static ret_t ret;
209
210 FILE *dirlist;
211
212 /* Check arguments... */
213 if (argc < 3)
214 ERRMSG("Give parameters: <ctl> <dirlist>");
215
216 /* Measure CPU-time... */
217 TIMER("total", 1);
218
219 /* Read control parameters... */
220 read_ctl(argc, argv, &ctl);
221 read_ret(argc, argv, &ctl, &ret);
222
223 /* Open directory list... */
224 if (!(dirlist = fopen(argv[2], "r")))
225 ERRMSG("Cannot open directory list!");
226
227 /* Loop over directories... */
228 while (fscanf(dirlist, "%4999s", ret.dir) != EOF) {
229
230 /* Write info... */
231 LOG(1, "\nRetrieve in directory %s...\n", ret.dir);
232
233 /* Read atmospheric data... */
234 read_atm(ret.dir, "atm_apr.tab", &ctl, &atm_apr);
235
236 /* Read observation data... */
237 read_obs(ret.dir, "obs_meas.tab", &ctl, &obs_meas);
238
239 /* Run retrieval... */
240 optimal_estimation(&ret, &ctl, &obs_meas, &obs_i, &atm_apr, &atm_i);
241
242 /* Measure CPU-time... */
243 TIMER("total", 2);
244 }
245
246 /* Write info... */
247 LOG(1, "\nRetrieval done...");
248
249 /* Measure CPU-time... */
250 TIMER("total", 3);
251
252 return EXIT_SUCCESS;
253}
void read_ctl(int argc, char *argv[], ctl_t *ctl)
Read forward model control parameters.
Definition: jurassic.c:4547
void read_atm(const char *dirname, const char *filename, const ctl_t *ctl, atm_t *atm)
Read atmospheric data.
Definition: jurassic.c:4442
void read_obs(const char *dirname, const char *filename, const ctl_t *ctl, obs_t *obs)
Read observation data.
Definition: jurassic.c:4697
#define TIMER(name, mode)
Start or stop a timer.
Definition: jurassic.h:201
void optimal_estimation(ret_t *ret, ctl_t *ctl, obs_t *obs_meas, obs_t *obs_i, atm_t *atm_apr, atm_t *atm_i)
Carry out optimal estimation retrieval.
Definition: retrieval.c:448
void read_ret(int argc, char *argv[], ctl_t *ctl, ret_t *ret)
Read retrieval control parameters.
Definition: retrieval.c:753
Forward model control parameters.
Definition: jurassic.h:541
Retrieval control parameters.
Definition: retrieval.c:32
Here is the call graph for this function: