37 int nc_result=(cmd); \
38 if(nc_result!=NC_NOERR) \
39 ERRMSG("%s", nc_strerror(nc_result)); \
134 double err_formod[ND];
137 double err_noise[ND];
186 const char *longname,
196 gsl_matrix * s_a_inv,
197 gsl_vector * sig_eps_inv);
261 gsl_vector * sig_noise,
262 gsl_vector * sig_formod,
263 gsl_vector * sig_eps_inv);
274 static atm_t atm_apr, atm_clim, atm_i;
275 static obs_t obs_i, obs_meas;
279 static FILE *in, *out;
281 static char filename[LEN], filename2[2 * LEN];
286 static int channel[ND], n, ntask = -1, rank, size;
293 MPI_Init(&argc, &argv);
294 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
295 MPI_Comm_size(MPI_COMM_WORLD, &size);
302 ERRMSG(
"Give parameters: <ctl> <filelist>");
305 read_ctl(argc, argv, &ctl);
309 tbl_t *tbl = read_tbl(&ctl);
312 const int nz = (int) scan_ctl(argc, argv,
"NZ", -1,
"", NULL);
314 ERRMSG(
"Too many altitudes!");
315 for (
int iz = 0; iz < nz; iz++)
316 z[iz] = scan_ctl(argc, argv,
"Z", iz,
"", NULL);
319 const int track0 = (int) scan_ctl(argc, argv,
"TRACK_MIN", -1,
"0", NULL);
320 const int track1 = (int) scan_ctl(argc, argv,
"TRACK_MAX", -1,
"134", NULL);
323 const int xtrack0 = (int) scan_ctl(argc, argv,
"XTRACK_MIN", -1,
"0", NULL);
325 (int) scan_ctl(argc, argv,
"XTRACK_MAX", -1,
"89", NULL);
328 const double sx = scan_ctl(argc, argv,
"SX", -1,
"8", NULL);
329 const double sy = scan_ctl(argc, argv,
"SY", -1,
"2", NULL);
336 printf(
"Read filelist: %s\n", argv[2]);
337 if (!(in = fopen(argv[2],
"r")))
338 ERRMSG(
"Cannot open filelist!");
341 while (fscanf(in,
"%s", filename) != EOF) {
344 if ((++ntask) % size != rank)
348 printf(
"Retrieve file %s on rank %d of %d (with %d threads)...\n",
349 filename, rank + 1, size, omp_get_max_threads());
359 for (
int id = 0;
id < ctl.nd;
id++) {
362 if (fabs(ctl.nu[
id] - ncd.
l1_nu[i]) < 0.1)
365 ERRMSG(
"Cannot identify radiance channel!");
374 for (
int iz = 0; iz < nz; iz++)
375 atm_clim.z[iz] = z[iz];
376 climatology(&ctl, &atm_clim);
383 for (
int track = track0; track <= track1; track++) {
386 TIMER(
"retrieval", 1);
389 for (
int xtrack = xtrack0; xtrack <= xtrack1; xtrack++) {
393 obs_meas.time[0] = ncd.
l1_time[track][xtrack];
394 obs_meas.obsz[0] = ncd.
l1_sat_z[track];
397 obs_meas.vplon[0] = ncd.
l1_lon[track][xtrack];
398 obs_meas.vplat[0] = ncd.
l1_lat[track][xtrack];
399 for (
int id = 0;
id < ctl.nd;
id++)
400 obs_meas.rad[
id][0] = ncd.
l1_rad[track][xtrack][channel[
id]];
403 for (
int id = 0;
id < ctl.nd;
id++)
404 if (ctl.nu[
id] >= 2000)
405 obs_meas.rad[id][0] = GSL_NAN;
408 copy_atm(&ctl, &atm_apr, &atm_clim, 0);
409 for (
int ip = 0; ip < atm_apr.np; ip++) {
410 atm_apr.time[ip] = obs_meas.time[0];
411 atm_apr.lon[ip] = obs_meas.vplon[0];
412 atm_apr.lat[ip] = obs_meas.vplat[0];
416 init_l2(&ncd, track, xtrack, &ctl, &atm_apr);
420 &atm_apr, &atm_i, &chisq[track][xtrack]);
423 for (
int id = 0;
id < ctl.nd;
id++)
424 obs_meas.rad[
id][0] = obs_i.rad[
id][0]
425 = ncd.
l1_rad[track][xtrack][channel[
id]];
426 formod(&ctl, tbl, &atm_i, &obs_i);
430 ni[track][xtrack] = 0;
431 for (
int id = 0;
id < ctl.nd;
id++)
432 if (ctl.nu[
id] >= 2000 && gsl_finite(obs_meas.rad[
id][0])) {
434 (BRIGHT(obs_meas.rad[
id][0], ncd.
l1_nu[channel[
id]])
435 - BRIGHT(obs_i.rad[
id][0], ncd.
l1_nu[channel[
id]]));
438 ni[track][xtrack] /= n;
442 TIMER(
"retrieval", 3);
450 sprintf(filename2,
"%s.nlte.tab", filename);
453 printf(
"Write non-LTE data: %s\n", filename2);
454 if (!(out = fopen(filename2,
"w")))
455 ERRMSG(
"Cannot create file!");
459 "# $1 = time (seconds since 2000-01-01, 00:00 UTC)\n"
460 "# $2 = longitude [deg]\n"
461 "# $3 = latitude [deg]\n"
462 "# $4 = solar zenith angle [deg]\n"
463 "# $5 = non-LTE index [K]\n" "# $6 = chi^2 of retrieval fit\n");
466 for (
int track = track0; track <= track1; track++) {
468 for (
int xtrack = xtrack0; xtrack <= xtrack1; xtrack++)
469 fprintf(out,
"%.2f %g %g %g %g %g\n",
471 ncd.
l1_lon[track][xtrack],
472 ncd.
l1_lat[track][xtrack],
474 ncd.
l1_lat[track][xtrack]), ni[track][xtrack],
475 chisq[track][xtrack]);
482 printf(
"Retrieval finished on rank %d of %d!\n", rank, size);
495 printf(
"MEMORY_ATM = %g MByte\n", 4. *
sizeof(atm_t) / 1024. / 1024.);
496 printf(
"MEMORY_CTL = %g MByte\n", 1. *
sizeof(ctl_t) / 1024. / 1024.);
497 printf(
"MEMORY_NCD = %g MByte\n", 1. *
sizeof(
ncd_t) / 1024. / 1024.);
498 printf(
"MEMORY_OBS = %g MByte\n", 3. *
sizeof(atm_t) / 1024. / 1024.);
499 printf(
"MEMORY_RET = %g MByte\n", 1. *
sizeof(
ret_t) / 1024. / 1024.);
500 printf(
"MEMORY_TBL = %g MByte\n", 1. *
sizeof(tbl_t) / 1024. / 1024.);
503 printf(
"SIZE_TASKS = %d\n", size);
504 printf(
"SIZE_THREADS = %d\n", omp_get_max_threads());
518 const char *longname,
525 if (nc_inq_varid(ncid, varname, varid) != NC_NOERR) {
528 NC(nc_def_var(ncid, varname, type, ndims, dimid, varid));
532 (ncid, *varid,
"long_name", strlen(longname), longname));
535 NC(nc_put_att_text(ncid, *varid,
"units", strlen(unit), unit));
545 gsl_vector *sig_eps_inv) {
547 double chisq_a, chisq_m = 0;
550 const size_t m = dy->size;
551 const size_t n = dx->size;
554 gsl_vector *x_aux = gsl_vector_alloc(n);
555 gsl_vector *y_aux = gsl_vector_alloc(m);
559 for (
size_t i = 0; i < m; i++)
561 gsl_pow_2(gsl_vector_get(dy, i) * gsl_vector_get(sig_eps_inv, i));
562 gsl_blas_dgemv(CblasNoTrans, 1.0, s_a_inv, dx, 0.0, x_aux);
563 gsl_blas_ddot(dx, x_aux, &chisq_a);
566 gsl_vector_free(x_aux);
567 gsl_vector_free(y_aux);
570 return (chisq_m + chisq_a) / (double) m;
583 for (
int lay = 0; lay <
L2_NLAY; lay++) {
586 for (
int track = 0; track <
L2_NTRACK; track++)
587 for (
int xtrack = 0; xtrack <
L2_NXTRACK; xtrack++) {
590 help[track][xtrack] = 0;
594 for (
int track2 = 0; track2 <
L2_NTRACK; track2++)
595 for (
int xtrack2 = 0; xtrack2 <
L2_NXTRACK; xtrack2++)
596 if (gsl_finite(x[track2][xtrack2][lay])
597 && x[track2][xtrack2][lay] > 0) {
598 const double w = exp(-gsl_pow_2((xtrack - xtrack2) / cx)
599 - gsl_pow_2((track - track2) / cy));
600 help[track][xtrack] += w * x[track2][xtrack2][lay];
606 help[track][xtrack] /= wsum;
608 help[track][xtrack] = GSL_NAN;
612 for (
int track = 0; track <
L2_NTRACK; track++)
613 for (
int xtrack = 0; xtrack <
L2_NXTRACK; xtrack++)
614 x[track][xtrack][lay] = help[track][xtrack];
627 static atm_t atm_airs;
629 double k[NW], p, q[NG], t, zmax = 0, zmin = 1000;
637 for (
int lay = 0; lay <
L2_NLAY; lay++)
638 if (gsl_finite(ncd->
l2_z[track][xtrack][lay])) {
639 atm_airs.z[atm_airs.np] = ncd->
l2_z[track][xtrack][lay];
640 atm_airs.p[atm_airs.np] = ncd->
l2_p[lay];
641 atm_airs.t[atm_airs.np] = ncd->
l2_t[track][xtrack][lay];
642 if ((++atm_airs.np) > NP)
643 ERRMSG(
"Too many layers!");
647 if (atm_airs.np <= 0)
651 for (
int ip = 0; ip < atm_airs.np; ip++) {
652 zmax = GSL_MAX(zmax, atm_airs.z[ip]);
653 zmin = GSL_MIN(zmin, atm_airs.z[ip]);
657 for (
int ip = 0; ip < atm->np; ip++) {
660 intpol_atm(ctl, &atm_airs, atm->z[ip], &p, &t, q, k);
664 if (atm->z[ip] > zmax)
665 w = GSL_MAX(1 - (atm->z[ip] - zmax) / 50, 0);
666 if (atm->z[ip] < zmin)
667 w = GSL_MAX(1 - (zmin - atm->z[ip]) / 50, 0);
670 atm->t[ip] = w * t + (1 - w) * atm->t[ip];
671 atm->p[ip] = w * p + (1 - w) * atm->p[ip];
683 const size_t n = a->size1;
686 for (
size_t i = 0; i < n && diag; i++)
687 for (
size_t j = i + 1; j < n; j++)
688 if (gsl_matrix_get(a, i, j) != 0) {
695 for (
size_t i = 0; i < n; i++)
696 gsl_matrix_set(a, i, i, 1 / gsl_matrix_get(a, i, i));
700 gsl_linalg_cholesky_decomp(a);
701 gsl_linalg_cholesky_invert(a);
714 const size_t m = a->size1;
715 const size_t n = a->size2;
718 gsl_matrix *aux = gsl_matrix_alloc(m, n);
721 if (transpose == 1) {
724 for (
size_t i = 0; i < m; i++)
725 for (
size_t j = 0; j < n; j++)
726 gsl_matrix_set(aux, i, j,
727 gsl_vector_get(b, i) * gsl_matrix_get(a, i, j));
730 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, aux, aux, 0.0, c);
734 else if (transpose == 2) {
737 for (
size_t i = 0; i < m; i++)
738 for (
size_t j = 0; j < n; j++)
739 gsl_matrix_set(aux, i, j,
740 gsl_matrix_get(a, i, j) * gsl_vector_get(b, j));
743 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, aux, aux, 0.0, c);
747 gsl_matrix_free(aux);
762 static int ipa[N], iqa[N];
764 double disq = 0, lmpar = 0.001;
771 const size_t m = obs2y(ctl, obs_meas, NULL, NULL, NULL);
772 const size_t n = atm2x(ctl, atm_apr, NULL, iqa, ipa);
773 if (m == 0 || n == 0) {
779 gsl_matrix *a = gsl_matrix_alloc(n, n);
780 gsl_matrix *cov = gsl_matrix_alloc(n, n);
781 gsl_matrix *k_i = gsl_matrix_alloc(m, n);
782 gsl_matrix *s_a_inv = gsl_matrix_alloc(n, n);
784 gsl_vector *b = gsl_vector_alloc(n);
785 gsl_vector *dx = gsl_vector_alloc(n);
786 gsl_vector *dy = gsl_vector_alloc(m);
787 gsl_vector *sig_eps_inv = gsl_vector_alloc(m);
788 gsl_vector *sig_formod = gsl_vector_alloc(m);
789 gsl_vector *sig_noise = gsl_vector_alloc(m);
790 gsl_vector *x_a = gsl_vector_alloc(n);
791 gsl_vector *x_i = gsl_vector_alloc(n);
792 gsl_vector *x_step = gsl_vector_alloc(n);
793 gsl_vector *y_aux = gsl_vector_alloc(m);
794 gsl_vector *y_i = gsl_vector_alloc(m);
795 gsl_vector *y_m = gsl_vector_alloc(m);
798 copy_atm(ctl, atm_i, atm_apr, 0);
799 copy_obs(ctl, obs_i, obs_meas, 0);
800 formod(ctl, tbl, atm_i, obs_i);
803 atm2x(ctl, atm_apr, x_a, NULL, NULL);
804 atm2x(ctl, atm_i, x_i, NULL, NULL);
805 obs2y(ctl, obs_meas, y_m, NULL, NULL);
806 obs2y(ctl, obs_i, y_i, NULL, NULL);
813 set_cov_meas(ret, ctl, obs_meas, sig_noise, sig_formod, sig_eps_inv);
816 gsl_vector_memcpy(dx, x_i);
817 gsl_vector_sub(dx, x_a);
818 gsl_vector_memcpy(dy, y_m);
819 gsl_vector_sub(dy, y_i);
825 kernel(ctl, tbl, atm_i, obs_i, k_i);
832 for (
int it = 1; it <= ret->
conv_itmax; it++) {
835 double chisq_old = *chisq;
839 kernel(ctl, tbl, atm_i, obs_i, k_i);
846 for (
size_t i = 0; i < m; i++)
847 gsl_vector_set(y_aux, i, gsl_vector_get(dy, i)
848 * gsl_pow_2(gsl_vector_get(sig_eps_inv, i)));
849 gsl_blas_dgemv(CblasTrans, 1.0, k_i, y_aux, 0.0, b);
850 gsl_blas_dgemv(CblasNoTrans, -1.0, s_a_inv, dx, 1.0, b);
853 for (
int it2 = 0; it2 < 20; it2++) {
856 gsl_matrix_memcpy(a, s_a_inv);
857 gsl_matrix_scale(a, 1 + lmpar);
858 gsl_matrix_add(a, cov);
861 gsl_linalg_cholesky_decomp(a);
862 gsl_linalg_cholesky_solve(a, b, x_step);
865 gsl_vector_add(x_i, x_step);
866 copy_atm(ctl, atm_i, atm_apr, 0);
867 copy_obs(ctl, obs_i, obs_meas, 0);
868 x2atm(ctl, x_i, atm_i);
871 for (
int ip = 0; ip < atm_i->np; ip++) {
872 atm_i->p[ip] = GSL_MIN(GSL_MAX(atm_i->p[ip], 5e-7), 5e4);
873 atm_i->t[ip] = GSL_MIN(GSL_MAX(atm_i->t[ip], 100), 400);
874 for (
int ig = 0; ig < ctl->ng; ig++)
875 atm_i->q[ig][ip] = GSL_MIN(GSL_MAX(atm_i->q[ig][ip], 0), 1);
876 for (
int iw = 0; iw < ctl->nw; iw++)
877 atm_i->k[iw][ip] = GSL_MAX(atm_i->k[iw][ip], 0);
881 formod(ctl, tbl, atm_i, obs_i);
882 obs2y(ctl, obs_i, y_i, NULL, NULL);
885 gsl_vector_memcpy(dx, x_i);
886 gsl_vector_sub(dx, x_a);
887 gsl_vector_memcpy(dy, y_m);
888 gsl_vector_sub(dy, y_i);
894 if (*chisq > chisq_old) {
896 gsl_vector_sub(x_i, x_step);
904 gsl_blas_ddot(x_step, b, &disq);
917 gsl_matrix_free(cov);
918 gsl_matrix_free(k_i);
919 gsl_matrix_free(s_a_inv);
924 gsl_vector_free(sig_eps_inv);
925 gsl_vector_free(sig_formod);
926 gsl_vector_free(sig_noise);
927 gsl_vector_free(x_a);
928 gsl_vector_free(x_i);
929 gsl_vector_free(x_step);
930 gsl_vector_free(y_aux);
931 gsl_vector_free(y_i);
932 gsl_vector_free(y_m);
944 printf(
"Read netCDF file: %s\n", filename);
945 NC(nc_open(filename, NC_WRITE, &ncd->
ncid));
948 NC(nc_inq_varid(ncd->
ncid,
"l1_time", &varid));
950 NC(nc_inq_varid(ncd->
ncid,
"l1_lon", &varid));
951 NC(nc_get_var_double(ncd->
ncid, varid, ncd->
l1_lon[0]));
952 NC(nc_inq_varid(ncd->
ncid,
"l1_lat", &varid));
953 NC(nc_get_var_double(ncd->
ncid, varid, ncd->
l1_lat[0]));
954 NC(nc_inq_varid(ncd->
ncid,
"l1_sat_z", &varid));
956 NC(nc_inq_varid(ncd->
ncid,
"l1_sat_lon", &varid));
958 NC(nc_inq_varid(ncd->
ncid,
"l1_sat_lat", &varid));
960 NC(nc_inq_varid(ncd->
ncid,
"l1_nu", &varid));
961 NC(nc_get_var_double(ncd->
ncid, varid, ncd->
l1_nu));
962 NC(nc_inq_varid(ncd->
ncid,
"l1_rad", &varid));
963 NC(nc_get_var_float(ncd->
ncid, varid, ncd->
l1_rad[0][0]));
966 NC(nc_inq_varid(ncd->
ncid,
"l2_z", &varid));
967 NC(nc_get_var_double(ncd->
ncid, varid, ncd->
l2_z[0][0]));
968 NC(nc_inq_varid(ncd->
ncid,
"l2_press", &varid));
969 NC(nc_get_var_double(ncd->
ncid, varid, ncd->
l2_p));
970 NC(nc_inq_varid(ncd->
ncid,
"l2_temp", &varid));
971 NC(nc_get_var_double(ncd->
ncid, varid, ncd->
l2_t[0][0]));
984 (int) scan_ctl(argc, argv,
"KERNEL_RECOMP", -1,
"3", NULL);
985 ret->
conv_itmax = (int) scan_ctl(argc, argv,
"CONV_ITMAX", -1,
"30", NULL);
986 ret->
conv_dmin = scan_ctl(argc, argv,
"CONV_DMIN", -1,
"0.1", NULL);
988 for (
int id = 0;
id < ctl->nd;
id++)
989 ret->
err_formod[
id] = scan_ctl(argc, argv,
"ERR_FORMOD",
id,
"0", NULL);
991 for (
int id = 0;
id < ctl->nd;
id++)
992 ret->
err_noise[
id] = scan_ctl(argc, argv,
"ERR_NOISE",
id,
"0", NULL);
994 ret->
err_press = scan_ctl(argc, argv,
"ERR_PRESS", -1,
"0", NULL);
995 ret->
err_press_cz = scan_ctl(argc, argv,
"ERR_PRESS_CZ", -1,
"-999", NULL);
996 ret->
err_press_ch = scan_ctl(argc, argv,
"ERR_PRESS_CH", -1,
"-999", NULL);
998 ret->
err_temp = scan_ctl(argc, argv,
"ERR_TEMP", -1,
"0", NULL);
999 ret->
err_temp_cz = scan_ctl(argc, argv,
"ERR_TEMP_CZ", -1,
"-999", NULL);
1000 ret->
err_temp_ch = scan_ctl(argc, argv,
"ERR_TEMP_CH", -1,
"-999", NULL);
1002 for (
int ig = 0; ig < ctl->ng; ig++) {
1003 ret->
err_q[ig] = scan_ctl(argc, argv,
"ERR_Q", ig,
"0", NULL);
1004 ret->
err_q_cz[ig] = scan_ctl(argc, argv,
"ERR_Q_CZ", ig,
"-999", NULL);
1005 ret->
err_q_ch[ig] = scan_ctl(argc, argv,
"ERR_Q_CH", ig,
"-999", NULL);
1008 for (
int iw = 0; iw < ctl->nw; iw++) {
1009 ret->
err_k[iw] = scan_ctl(argc, argv,
"ERR_K", iw,
"0", NULL);
1010 ret->
err_k_cz[iw] = scan_ctl(argc, argv,
"ERR_K_CZ", iw,
"-999", NULL);
1011 ret->
err_k_ch[iw] = scan_ctl(argc, argv,
"ERR_K_CH", iw,
"-999", NULL);
1026 const size_t n = s_a->size1;
1029 gsl_vector *x_a = gsl_vector_alloc(n);
1032 atm2x(ctl, atm, x_a, NULL, NULL);
1033 for (
size_t i = 0; i < n; i++) {
1035 gsl_vector_set(x_a, i, ret->
err_press / 100 * gsl_vector_get(x_a, i));
1037 gsl_vector_set(x_a, i, ret->
err_temp);
1038 for (
int ig = 0; ig < ctl->ng; ig++)
1039 if (iqa[i] == IDXQ(ig))
1040 gsl_vector_set(x_a, i, ret->
err_q[ig] / 100 * gsl_vector_get(x_a, i));
1041 for (
int iw = 0; iw < ctl->nw; iw++)
1042 if (iqa[i] == IDXK(iw))
1043 gsl_vector_set(x_a, i, ret->
err_k[iw]);
1047 for (
size_t i = 0; i < n; i++)
1048 if (gsl_pow_2(gsl_vector_get(x_a, i)) <= 0)
1049 ERRMSG(
"Check a priori data (zero standard deviation)!");
1052 gsl_matrix_set_zero(s_a);
1053 for (
size_t i = 0; i < n; i++)
1054 gsl_matrix_set(s_a, i, i, gsl_pow_2(gsl_vector_get(x_a, i)));
1057 for (
size_t i = 0; i < n; i++)
1058 for (
size_t j = 0; j < n; j++)
1059 if (i != j && iqa[i] == iqa[j]) {
1062 double cz = 0, ch = 0;
1065 if (iqa[i] == IDXP) {
1071 if (iqa[i] == IDXT) {
1077 for (
int ig = 0; ig < ctl->ng; ig++)
1078 if (iqa[i] == IDXQ(ig)) {
1084 for (
int iw = 0; iw < ctl->nw; iw++)
1085 if (iqa[i] == IDXK(iw)) {
1091 if (cz > 0 && ch > 0) {
1094 double x0[3], x1[3];
1095 geo2cart(0, atm->lon[ipa[i]], atm->lat[ipa[i]], x0);
1096 geo2cart(0, atm->lon[ipa[j]], atm->lat[ipa[j]], x1);
1100 exp(-DIST(x0, x1) / ch -
1101 fabs(atm->z[ipa[i]] - atm->z[ipa[j]]) / cz);
1104 gsl_matrix_set(s_a, i, j, gsl_vector_get(x_a, i)
1105 * gsl_vector_get(x_a, j) * rho);
1110 gsl_vector_free(x_a);
1119 gsl_vector *sig_noise,
1120 gsl_vector *sig_formod,
1121 gsl_vector *sig_eps_inv) {
1123 static obs_t obs_err;
1126 const size_t m = sig_eps_inv->size;
1129 copy_obs(ctl, &obs_err, obs, 1);
1130 for (
int ir = 0; ir < obs_err.nr; ir++)
1131 for (
int id = 0;
id < ctl->nd;
id++)
1133 = (gsl_finite(obs->rad[
id][ir]) ? ret->
err_noise[
id] : GSL_NAN);
1134 obs2y(ctl, &obs_err, sig_noise, NULL, NULL);
1137 copy_obs(ctl, &obs_err, obs, 1);
1138 for (
int ir = 0; ir < obs_err.nr; ir++)
1139 for (
int id = 0;
id < ctl->nd;
id++)
1141 = fabs(ret->
err_formod[
id] / 100 * obs->rad[
id][ir]);
1142 obs2y(ctl, &obs_err, sig_formod, NULL, NULL);
1145 for (
size_t i = 0; i < m; i++)
1146 gsl_vector_set(sig_eps_inv, i,
1147 1 / sqrt(gsl_pow_2(gsl_vector_get(sig_noise, i))
1148 + gsl_pow_2(gsl_vector_get(sig_formod, i))));
1151 for (
size_t i = 0; i < m; i++)
1152 if (gsl_vector_get(sig_eps_inv, i) <= 0)
1153 ERRMSG(
"Check measurement errors (zero standard deviation)!");
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[])
#define NC(cmd)
Execute netCDF library command and check result.
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 fill_gaps(double x[L2_NTRACK][L2_NXTRACK][L2_NLAY], double cx, double cy)
Fill data gaps in L2 data.
#define L1_NXTRACK
Across-track size of AIRS radiance granule (don't change).
#define L1_NTRACK
Along-track size of AIRS radiance granule (don't change).
#define L2_NXTRACK
Across-track size of AIRS retrieval granule (don't change).
void read_nc(char *filename, ncd_t *ncd)
Read netCDF file.
void add_var(int ncid, const char *varname, const char *unit, const char *longname, int type, int dimid[], int *varid, int ndims)
Create variable in netCDF file.
#define L2_NLAY
Number of AIRS pressure layers (don't change).
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, double *chisq)
Carry out optimal estimation retrieval.
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 read_ret_ctl(int argc, char *argv[], ctl_t *ctl, ret_t *ret)
Read retrieval control parameters.
void matrix_invert(gsl_matrix *a)
Invert symmetric matrix.
#define L1_NCHAN
Number of AIRS radiance channels (don't change).
#define L2_NTRACK
Along-track size of AIRS retrieval granule (don't change).
double cost_function(gsl_vector *dx, gsl_vector *dy, gsl_matrix *s_a_inv, gsl_vector *sig_eps_inv)
Compute cost function.
void init_l2(ncd_t *ncd, int track, int xtrack, ctl_t *ctl, atm_t *atm)
Initialize with AIRS Level-2 data.
double l1_lat[L1_NTRACK][L1_NXTRACK]
Footprint latitude [deg].
double l1_lon[L1_NTRACK][L1_NXTRACK]
Footprint longitude [deg].
double l2_z[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Altitude [km].
double l1_sat_lat[L1_NTRACK]
Satellite latitude [deg].
double l1_sat_lon[L1_NTRACK]
Satellite longitude [deg].
double l2_p[L2_NLAY]
Pressure [hPa].
double l2_t[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Temperature [K].
double l1_nu[L1_NCHAN]
Channel frequencies [cm^-1].
double l1_sat_z[L1_NTRACK]
Satellite altitude [km].
float l1_rad[L1_NTRACK][L1_NXTRACK][L1_NCHAN]
Radiance [W/(m^2 sr cm^-1)].
double l1_time[L1_NTRACK][L1_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
double err_k_ch[NW]
Horizontal correlation length for extinction error [km].
double err_press_cz
Vertical correlation length for pressure error [km].
double err_k_cz[NW]
Vertical correlation length for extinction error [km].
double err_press
Pressure error [%].
double err_q_ch[NG]
Horizontal correlation length for volume mixing ratio error [km].
double err_noise[ND]
Noise error [W/(m^2 sr cm^-1)].
double err_formod[ND]
Forward model error [%].
double err_q_cz[NG]
Vertical correlation length for volume mixing ratio error [km].
double err_temp_cz
Vertical correlation length for temperature error [km].
double conv_dmin
Minimum normalized step size in state space.
double err_temp
Temperature error [K].
double err_temp_ch
Horizontal correlation length for temperature error [km].
int kernel_recomp
Recomputation of kernel matrix (number of iterations).
int conv_itmax
Maximum number of iterations.
double err_press_ch
Horizontal correlation length for pressure error [km].
double err_q[NG]
Volume mixing ratio error [%].
double err_k[NW]
Extinction error [1/km].