50 double err_formod[
ND];
141 gsl_matrix * s_a_inv,
142 gsl_vector * sig_eps_inv);
185 gsl_vector * sig_noise,
186 gsl_vector * sig_formod,
187 gsl_vector * sig_eps_inv);
191 const char *quantity,
205 static atm_t atm_i, atm_apr;
207 static obs_t obs_i, obs_meas;
214 ERRMSG(
"Give parameters: <ctl> <dirlist>");
224 if (!(dirlist = fopen(argv[2],
"r")))
225 ERRMSG(
"Cannot open directory list!");
228 while (fscanf(dirlist,
"%4999s", ret.
dir) != EOF) {
231 LOG(1,
"\nRetrieve in directory %s...\n", ret.
dir);
237 read_obs(ret.
dir,
"obs_meas.tab", &ctl, &obs_meas);
247 LOG(1,
"\nRetrieval done...");
265 static atm_t atm_cont, atm_res;
267 size_t i, n0[
NQ], n1[
NQ];
270 const size_t n = avk->size1;
273 for (
int iq = 0; iq <
NQ; iq++) {
275 for (i = 0; i < n; i++) {
276 if (iqa[i] == iq && n0[iq] ==
N)
279 n1[iq] = i - n0[iq] + 1;
290 for (
int ig = 0; ig < ctl->
ng; ig++)
292 atm_cont.
q[ig], atm_res.
q[ig]);
293 for (
int iw = 0; iw < ctl->
nw; iw++)
295 atm_cont.
k[iw], atm_res.
k[iw]);
299 for (
int icl = 0; icl < ctl->
ncl; icl++)
301 &atm_cont.
clk[icl], &atm_res.
clk[icl]);
305 for (
int isf = 0; isf < ctl->
nsf; isf++)
327 for (
size_t i = 0; i < n1[iq]; i++) {
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);
334 res[ipa[n0[iq] + i]] = 1 / gsl_matrix_get(avk, n0[iq] + i, n0[iq] + i);
344 gsl_vector *sig_eps_inv) {
346 double chisq_a, chisq_m = 0;
349 const size_t m = dy->size;
350 const size_t n = dx->size;
353 gsl_vector *x_aux = gsl_vector_alloc(n);
354 gsl_vector *y_aux = gsl_vector_alloc(m);
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);
364 gsl_vector_free(x_aux);
365 gsl_vector_free(y_aux);
368 return (chisq_m + chisq_a) / (double) m;
379 const size_t n = a->size1;
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) {
391 for (
size_t i = 0; i < n; i++)
392 gsl_matrix_set(a, i, i, 1 / gsl_matrix_get(a, i, i));
396 gsl_linalg_cholesky_decomp(a);
397 gsl_linalg_cholesky_invert(a);
410 const size_t m = a->size1;
411 const size_t n = a->size2;
414 gsl_matrix *aux = gsl_matrix_alloc(m, n);
417 if (transpose == 1) {
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));
426 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, aux, aux, 0.0, c);
430 else if (transpose == 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));
439 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, aux, aux, 0.0, c);
443 gsl_matrix_free(aux);
456 static int ipa[
N], iqa[
N];
460 char filename[2 *
LEN];
462 double chisq, disq = 0, lmpar = 0.001;
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!");
477 LOG(1,
"Problem size: m= %d / n= %d "
478 "(alloc= %.4g MB / stat= %.4g MB)",
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.);
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);
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);
507 formod(ctl, atm_i, obs_i);
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);
518 atm_i, obs_i,
"x",
"x",
"r");
522 set_cov_meas(ret, ctl, obs_meas, sig_noise, sig_formod, sig_eps_inv);
525 sprintf(filename,
"%s/costs.tab", ret->
dir);
526 if (!(out = fopen(filename,
"w")))
527 ERRMSG(
"Cannot create cost function file!");
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");
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);
546 LOG(1,
"it= %d / chi^2/m= %g", it, chisq);
549 fprintf(out,
"%d %g %d %d\n", it, chisq, (
int) m, (
int) n);
552 kernel(ctl, atm_i, obs_i, k_i);
562 double chisq_old = chisq;
566 kernel(ctl, atm_i, obs_i, k_i);
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);
580 for (
int it2 = 0; it2 < 20; it2++) {
583 gsl_matrix_memcpy(a, s_a_inv);
584 gsl_matrix_scale(a, 1 + lmpar);
585 gsl_matrix_add(a, cov);
588 gsl_linalg_cholesky_decomp(a);
589 gsl_linalg_cholesky_solve(a, b, x_step);
592 gsl_vector_add(x_i, x_step);
595 x2atm(ctl, x_i, atm_i);
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);
608 for (
int icl = 0; icl < ctl->
ncl; icl++)
609 atm_i->
clk[icl] =
MAX(atm_i->
clk[icl], 0);
613 for (
int isf = 0; isf < ctl->
nsf; isf++)
617 formod(ctl, atm_i, obs_i);
618 obs2y(ctl, obs_i, y_i, NULL, NULL);
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);
630 if (chisq > chisq_old) {
632 gsl_vector_sub(x_i, x_step);
640 LOG(1,
"it= %d / chi^2/m= %g", it, chisq);
643 fprintf(out,
"%d %g %d %d\n", it, chisq, (
int) m, (
int) n);
646 gsl_blas_ddot(x_step, b, &disq);
661 atm_i, obs_i,
"y",
"x",
"r");
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);
678 gsl_matrix_add(cov, s_a_inv);
683 atm_i, obs_i,
"x",
"x",
"r");
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)));
693 atm_i, obs_i,
"x",
"x",
"r");
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);
703 atm_i, obs_i,
"x",
"y",
"c");
715 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, gain, k_i, 0.0, a);
717 atm_i, obs_i,
"x",
"x",
"r");
723 gsl_matrix_free(auxnm);
724 gsl_matrix_free(corr);
725 gsl_matrix_free(gain);
733 gsl_matrix_free(cov);
734 gsl_matrix_free(k_i);
735 gsl_matrix_free(s_a_inv);
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);
761 (int)
scan_ctl(argc, argv,
"KERNEL_RECOMP", -1,
"3", NULL);
768 for (
int id = 0;
id < ctl->
nd;
id++)
771 for (
int id = 0;
id < ctl->
nd;
id++)
782 for (
int ig = 0; ig < ctl->
ng; ig++) {
783 ret->
err_q[ig] =
scan_ctl(argc, argv,
"ERR_Q", ig,
"0", NULL);
788 for (
int iw = 0; iw < ctl->
nw; iw++) {
789 ret->
err_k[iw] =
scan_ctl(argc, argv,
"ERR_K", iw,
"0", NULL);
796 for (
int icl = 0; icl < ctl->
ncl; icl++)
802 for (
int isf = 0; isf < ctl->
nsf; isf++)
819 const size_t n = s_a->size1;
822 gsl_vector *x_a = gsl_vector_alloc(n);
825 atm2x(ctl, atm, x_a, NULL, NULL);
826 for (
size_t i = 0; i < n; i++) {
828 gsl_vector_set(x_a, i, ret->
err_press / 100 * gsl_vector_get(x_a, i));
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]);
838 gsl_vector_set(x_a, i, ret->
err_clz);
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]);
845 gsl_vector_set(x_a, i, ret->
err_sfz);
847 gsl_vector_set(x_a, i, ret->
err_sfp);
849 gsl_vector_set(x_a, i, ret->
err_sft);
850 for (
int isf = 0; isf < ctl->
nsf; isf++)
852 gsl_vector_set(x_a, i, ret->
err_sfeps[isf]);
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)!");
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)));
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]) {
875 if (iqa[i] ==
IDXP) {
881 if (iqa[i] ==
IDXT) {
887 for (
int ig = 0; ig < ctl->
ng; ig++)
888 if (iqa[i] ==
IDXQ(ig)) {
894 for (
int iw = 0; iw < ctl->
nw; iw++)
895 if (iqa[i] ==
IDXK(iw)) {
901 if (cz > 0 && ch > 0) {
909 exp(-
DIST(x0, x1) / ch -
910 fabs(atm->
z[ipa[i]] - atm->
z[ipa[j]]) / cz);
913 gsl_matrix_set(s_a, i, j, gsl_vector_get(x_a, i)
914 * gsl_vector_get(x_a, j) * rho);
919 gsl_vector_free(x_a);
928 gsl_vector *sig_noise,
929 gsl_vector *sig_formod,
930 gsl_vector *sig_eps_inv) {
932 static obs_t obs_err;
935 const size_t m = sig_eps_inv->size;
939 for (
int ir = 0; ir < obs_err.
nr; ir++)
940 for (
int id = 0;
id < ctl->
nd;
id++)
942 = (isfinite(obs->
rad[
id][ir]) ? ret->
err_noise[
id] : NAN);
943 obs2y(ctl, &obs_err, sig_noise, NULL, NULL);
947 for (
int ir = 0; ir < obs_err.
nr; ir++)
948 for (
int id = 0;
id < ctl->
nd;
id++)
951 obs2y(ctl, &obs_err, sig_formod, NULL, NULL);
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))
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)!");
969 const char *quantity,
975 static atm_t atm_aux;
980 const size_t n = s->size1;
983 gsl_vector *x_aux = gsl_vector_alloc(n);
986 for (
size_t i = 0; i < n; i++)
987 gsl_vector_set(x_aux, i, sqrt(gsl_matrix_get(s, i, i)));
991 x2atm(ctl, x_aux, &atm_aux);
992 sprintf(filename,
"atm_err_%s.tab", quantity);
996 gsl_vector_free(x_aux);
void write_atm(const char *dirname, const char *filename, const ctl_t *ctl, const atm_t *atm)
Write atmospheric data.
void read_ctl(int argc, char *argv[], ctl_t *ctl)
Read forward model control parameters.
void x2atm(const ctl_t *ctl, const gsl_vector *x, atm_t *atm)
Decompose parameter vector or state vector.
void formod(const ctl_t *ctl, atm_t *atm, obs_t *obs)
Determine ray paths and compute radiative transfer.
void read_atm(const char *dirname, const char *filename, const ctl_t *ctl, atm_t *atm)
Read atmospheric data.
void copy_obs(const ctl_t *ctl, obs_t *obs_dest, const obs_t *obs_src, const int init)
Copy and initialize observation data.
void read_obs(const char *dirname, const char *filename, const ctl_t *ctl, obs_t *obs)
Read observation data.
void write_obs(const char *dirname, const char *filename, const ctl_t *ctl, const obs_t *obs)
Write observation data.
double scan_ctl(int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
void copy_atm(const ctl_t *ctl, atm_t *atm_dest, const atm_t *atm_src, const int init)
Copy and initialize atmospheric data.
size_t obs2y(const ctl_t *ctl, const obs_t *obs, gsl_vector *y, int *ida, int *ira)
Compose measurement vector.
size_t atm2x(const ctl_t *ctl, const atm_t *atm, gsl_vector *x, int *iqa, int *ipa)
Compose state vector or parameter vector.
void geo2cart(const double z, const double lon, const double lat, double *x)
Convert geolocation to Cartesian coordinates.
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.
void kernel(ctl_t *ctl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute Jacobians.
JURASSIC library declarations.
#define N
Maximum size of state vector.
#define IDXCLZ
Index for cloud layer height.
#define LEN
Maximum length of ASCII data lines.
#define POW2(x)
Compute x^2.
#define MIN(a, b)
Macro to determine the minimum of two values.
#define IDXCLDZ
Index for cloud layer depth.
#define ERRMSG(...)
Print error message and quit program.
#define IDXK(iw)
Indices for extinction.
#define NQ
Maximum number of quantities.
#define ND
Maximum number of radiance channels.
#define IDXSFT
Index for surface layer temperature.
#define IDXSFEPS(isf)
Indices for surface layer emissivity.
#define IDXCLK(icl)
Indices for cloud layer extinction.
#define TIMER(name, mode)
Start or stop a timer.
#define IDXP
Index for pressure.
#define IDXSFP
Index for surface layer pressure.
#define NG
Maximum number of emitters.
#define LOG(level,...)
Print log message.
#define NSF
Maximum number of surface layer spectral grid points.
#define NCL
Maximum number of cloud layer spectral grid points.
#define DIST(a, b)
Compute Cartesian distance between two vectors.
#define IDXQ(ig)
Indices for volume mixing ratios.
#define NW
Maximum number of spectral windows.
#define MAX(a, b)
Macro to determine the maximum of two values.
#define IDXT
Index for temperature.
#define IDXSFZ
Index for surface layer height.
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.
int main(int argc, char *argv[])
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.
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.
void analyze_avk(ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *avk)
Compute information content and resolution.
void read_ret(int argc, char *argv[], ctl_t *ctl, ret_t *ret)
Read retrieval control parameters.
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.
void write_stddev(const char *quantity, ret_t *ret, ctl_t *ctl, atm_t *atm, gsl_matrix *s)
Write retrieval error to file.
void matrix_invert(gsl_matrix *a)
Invert symmetric matrix.
double cost_function(gsl_vector *dx, gsl_vector *dy, gsl_matrix *s_a_inv, gsl_vector *sig_eps_inv)
Compute cost function.
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.
double sfeps[NSF]
Surface emissivity.
double k[NW][NP]
Extinction [km^-1].
double sfz
Surface height [km].
double lat[NP]
Latitude [deg].
double lon[NP]
Longitude [deg].
double t[NP]
Temperature [K].
double sfp
Surface pressure [hPa].
double clz
Cloud layer height [km].
int np
Number of data points.
double cldz
Cloud layer depth [km].
double sft
Surface temperature [K].
double z[NP]
Altitude [km].
double clk[NCL]
Cloud layer extinction [km^-1].
double q[NG][NP]
Volume mixing ratio [ppv].
double p[NP]
Pressure [hPa].
Forward model control parameters.
int nw
Number of spectral windows.
int ng
Number of emitters.
int nd
Number of radiance channels.
int ncl
Number of cloud layer spectral grid points.
int nsf
Number of surface layer spectral grid points.
Observation geometry and radiance data.
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
int nr
Number of ray paths.
Retrieval control parameters.
double err_sfz
Surface height error [km].
double err_press_cz
Vertical correlation length for pressure error [km].
double err_sfp
Surface pressure error [hPa].
double err_press
Pressure error [%].
double err_k_cz[NW]
Vertical correlation length for extinction error [km].
int err_ana
Carry out error analysis (0=no, 1=yes).
double err_k_ch[NW]
Horizontal correlation length for extinction error [km].
double err_temp_cz
Vertical correlation length for temperature error [km].
double err_formod[ND]
Forward model error [%].
double conv_dmin
Minimum normalized step size in state space.
double err_temp
Temperature error [K].
double err_q_cz[NG]
Vertical correlation length for volume mixing ratio error [km].
double err_noise[ND]
Noise error [W/(m^2 sr cm^-1)].
double err_clz
Cloud height error [km].
double err_sft
Surface temperature error [K].
double err_clk[NCL]
Cloud extinction error [km^-1].
double err_temp_ch
Horizontal correlation length for temperature error [km].
int kernel_recomp
Re-computation of kernel matrix (number of iterations).
int conv_itmax
Maximum number of iterations.
double err_cldz
Cloud depth error [km].
double err_press_ch
Horizontal correlation length for pressure error [km].
double err_q_ch[NG]
Horizontal correlation length for volume mixing ratio error [km].
double err_q[NG]
Volume mixing ratio error [%].
char dir[LEN]
Working directory.
double err_sfeps[NSF]
Surface emissivity error.
double err_k[NW]
Extinction error [km^-1].