Retrieval of non-LTE index. More...
#include <mpi.h>#include <omp.h>#include <netcdf.h>#include "jurassic.h"Go to the source code of this file.
Data Structures | |
| struct | ncd_t |
| Buffer for netCDF data. More... | |
Macros | |
| #define | NC(cmd) |
| Execute netCDF library command and check result. More... | |
| #define | L1_NCHAN 34 |
| Number of AIRS radiance channels (don't change). More... | |
| #define | L1_NTRACK 135 |
| Along-track size of AIRS radiance granule (don't change). More... | |
| #define | L1_NXTRACK 90 |
| Across-track size of AIRS radiance granule (don't change). More... | |
| #define | L2_NLAY 27 |
| Number of AIRS pressure layers (don't change). More... | |
| #define | L2_NTRACK 45 |
| Along-track size of AIRS retrieval granule (don't change). More... | |
| #define | L2_NXTRACK 30 |
| Across-track size of AIRS retrieval granule (don't change). More... | |
Functions | |
| 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. More... | |
| void | fill_gaps (double x[L2_NTRACK][L2_NXTRACK][L2_NLAY], double cx, double cy) |
| Fill data gaps in L2 data. More... | |
| void | init_l2 (ncd_t *ncd, int track, int xtrack, ctl_t *ctl, atm_t *atm) |
| Initialize with AIRS Level-2 data. More... | |
| 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. More... | |
| void | read_nc (char *filename, ncd_t *ncd) |
| Read netCDF file. More... | |
| int | main (int argc, char *argv[]) |
Detailed Description
Retrieval of non-LTE index.
Definition in file nlte.c.
Macro Definition Documentation
◆ NC
| #define NC | ( | cmd | ) |
◆ L1_NCHAN
| #define L1_NCHAN 34 |
◆ L1_NTRACK
| #define L1_NTRACK 135 |
◆ L1_NXTRACK
| #define L1_NXTRACK 90 |
◆ L2_NLAY
| #define L2_NLAY 27 |
◆ L2_NTRACK
| #define L2_NTRACK 45 |
◆ L2_NXTRACK
| #define L2_NXTRACK 30 |
Function Documentation
◆ add_var()
| 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.
Add variable to netCDF file.
Definition at line 415 of file nlte.c.
423 {
424
425 /* Check if variable exists... */
426 if (nc_inq_varid(ncid, varname, varid) != NC_NOERR) {
427
428 /* Define variable... */
429 NC(nc_def_var(ncid, varname, type, ndims, dimid, varid));
430
431 /* Set long name... */
432 NC(nc_put_att_text
433 (ncid, *varid, "long_name", strlen(longname), longname));
434
435 /* Set units... */
437 }
438}
◆ fill_gaps()
| void fill_gaps | ( | double | x[L2_NTRACK][L2_NXTRACK][L2_NLAY], |
| double | cx, | ||
| double | cy | ||
| ) |
Fill data gaps in L2 data.
Definition at line 442 of file nlte.c.
445 {
446
448
449 /* Loop over layers... */
451
452 /* Loop over grid points... */
455
456 /* Init... */
457 help[track][xtrack] = 0;
458 double wsum = 0;
459
460 /* Averrage data points... */
463 if (gsl_finite(x[track2][xtrack2][lay])
464 && x[track2][xtrack2][lay] > 0) {
465 const double w = exp(-gsl_pow_2((xtrack - xtrack2) / cx)
466 - gsl_pow_2((track - track2) / cy));
467 help[track][xtrack] += w * x[track2][xtrack2][lay];
468 wsum += w;
469 }
470
471 /* Normalize... */
472 if (wsum > 0)
473 help[track][xtrack] /= wsum;
474 else
475 help[track][xtrack] = GSL_NAN;
476 }
477
478 /* Copy grid points... */
481 x[track][xtrack][lay] = help[track][xtrack];
482 }
483}
◆ init_l2()
| void init_l2 | ( | ncd_t * | ncd, |
| int | track, | ||
| int | xtrack, | ||
| ctl_t * | ctl, | ||
| atm_t * | atm | ||
| ) |
Initialize with AIRS Level-2 data.
Definition at line 487 of file nlte.c.
492 {
493
494 static atm_t atm_airs;
495
496 double k[NW], p, q[NG], t, zmax = 0, zmin = 1000;
497
498 /* Reset track- and xtrack-index to match Level-2 data... */
499 track /= 3;
500 xtrack /= 3;
501
502 /* Store AIRS data in atmospheric data struct... */
503 atm_airs.np = 0;
506 atm_airs.z[atm_airs.np] = ncd->l2_z[track][xtrack][lay];
507 atm_airs.p[atm_airs.np] = ncd->l2_p[lay];
508 atm_airs.t[atm_airs.np] = ncd->l2_t[track][xtrack][lay];
509 if ((++atm_airs.np) > NP)
510 ERRMSG("Too many layers!");
511 }
512
513 /* Check number of levels... */
514 if (atm_airs.np <= 0)
515 return;
516
517 /* Get height range of AIRS data... */
518 for (int ip = 0; ip < atm_airs.np; ip++) {
519 zmax = GSL_MAX(zmax, atm_airs.z[ip]);
520 zmin = GSL_MIN(zmin, atm_airs.z[ip]);
521 }
522
523 /* Merge AIRS data... */
524 for (int ip = 0; ip < atm->np; ip++) {
525
526 /* Interpolate AIRS data... */
527 intpol_atm(ctl, &atm_airs, atm->z[ip], &p, &t, q, k);
528
529 /* Weighting factor... */
530 double w = 1;
531 if (atm->z[ip] > zmax)
532 w = GSL_MAX(1 - (atm->z[ip] - zmax) / 50, 0);
533 if (atm->z[ip] < zmin)
534 w = GSL_MAX(1 - (zmin - atm->z[ip]) / 50, 0);
535
536 /* Merge... */
537 atm->t[ip] = w * t + (1 - w) * atm->t[ip];
538 atm->p[ip] = w * p + (1 - w) * atm->p[ip];
539 }
540}
◆ 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, | ||
| double * | chisq | ||
| ) |
Carry out optimal estimation retrieval.
Definition at line 544 of file nlte.c.
552 {
553
554 static int ipa[N], iqa[N];
555
556 double disq = 0, lmpar = 0.001;
557
558 /* ------------------------------------------------------------
559 Initialize...
560 ------------------------------------------------------------ */
561
562 /* Get sizes... */
563 const size_t m = obs2y(ctl, obs_meas, NULL, NULL, NULL);
564 const size_t n = atm2x(ctl, atm_apr, NULL, iqa, ipa);
565 if (m == 0 || n == 0) {
566 *chisq = GSL_NAN;
567 return;
568 }
569
570 /* Allocate... */
571 gsl_matrix *a = gsl_matrix_alloc(n, n);
572 gsl_matrix *cov = gsl_matrix_alloc(n, n);
573 gsl_matrix *k_i = gsl_matrix_alloc(m, n);
574 gsl_matrix *s_a_inv = gsl_matrix_alloc(n, n);
575
576 gsl_vector *b = gsl_vector_alloc(n);
577 gsl_vector *dx = gsl_vector_alloc(n);
578 gsl_vector *dy = gsl_vector_alloc(m);
579 gsl_vector *sig_eps_inv = gsl_vector_alloc(m);
580 gsl_vector *sig_formod = gsl_vector_alloc(m);
581 gsl_vector *sig_noise = gsl_vector_alloc(m);
582 gsl_vector *x_a = gsl_vector_alloc(n);
583 gsl_vector *x_i = gsl_vector_alloc(n);
584 gsl_vector *x_step = gsl_vector_alloc(n);
585 gsl_vector *y_aux = gsl_vector_alloc(m);
586 gsl_vector *y_i = gsl_vector_alloc(m);
587 gsl_vector *y_m = gsl_vector_alloc(m);
588
589 /* Set initial state... */
590 copy_atm(ctl, atm_i, atm_apr, 0);
591 copy_obs(ctl, obs_i, obs_meas, 0);
592 formod(ctl, tbl, atm_i, obs_i);
593
594 /* Set state vectors and observation vectors... */
595 atm2x(ctl, atm_apr, x_a, NULL, NULL);
596 atm2x(ctl, atm_i, x_i, NULL, NULL);
597 obs2y(ctl, obs_meas, y_m, NULL, NULL);
598 obs2y(ctl, obs_i, y_i, NULL, NULL);
599
600 /* Set inverse a priori covariance S_a^-1... */
601 set_cov_apr(ret, ctl, atm_apr, iqa, ipa, s_a_inv);
602 matrix_invert(s_a_inv);
603
604 /* Get measurement errors... */
605 set_cov_meas(ret, ctl, obs_meas, sig_noise, sig_formod, sig_eps_inv);
606
607 /* Determine dx = x_i - x_a and dy = y - F(x_i) ... */
608 gsl_vector_memcpy(dx, x_i);
609 gsl_vector_sub(dx, x_a);
610 gsl_vector_memcpy(dy, y_m);
611 gsl_vector_sub(dy, y_i);
612
613 /* Compute cost function... */
614 *chisq = cost_function(dx, dy, s_a_inv, sig_eps_inv);
615
616 /* Compute initial kernel... */
617 kernel(ctl, tbl, atm_i, obs_i, k_i);
618
619 /* ------------------------------------------------------------
620 Levenberg-Marquardt minimization...
621 ------------------------------------------------------------ */
622
623 /* Outer loop... */
624 for (int it = 1; it <= ret->conv_itmax; it++) {
625
626 /* Store current cost function value... */
627 double chisq_old = *chisq;
628
629 /* Compute kernel matrix K_i... */
630 if (it > 1 && it % ret->kernel_recomp == 0)
631 kernel(ctl, tbl, atm_i, obs_i, k_i);
632
633 /* Compute K_i^T * S_eps^{-1} * K_i ... */
634 if (it == 1 || it % ret->kernel_recomp == 0)
635 matrix_product(k_i, sig_eps_inv, 1, cov);
636
637 /* Determine b = K_i^T * S_eps^{-1} * dy - S_a^{-1} * dx ... */
638 for (size_t i = 0; i < m; i++)
639 gsl_vector_set(y_aux, i, gsl_vector_get(dy, i)
640 * gsl_pow_2(gsl_vector_get(sig_eps_inv, i)));
641 gsl_blas_dgemv(CblasTrans, 1.0, k_i, y_aux, 0.0, b);
642 gsl_blas_dgemv(CblasNoTrans, -1.0, s_a_inv, dx, 1.0, b);
643
644 /* Inner loop... */
645 for (int it2 = 0; it2 < 20; it2++) {
646
647 /* Compute A = (1 + lmpar) * S_a^{-1} + K_i^T * S_eps^{-1} * K_i ... */
648 gsl_matrix_memcpy(a, s_a_inv);
649 gsl_matrix_scale(a, 1 + lmpar);
650 gsl_matrix_add(a, cov);
651
652 /* Solve A * x_step = b by means of Cholesky decomposition... */
653 gsl_linalg_cholesky_decomp(a);
654 gsl_linalg_cholesky_solve(a, b, x_step);
655
656 /* Update atmospheric state... */
657 gsl_vector_add(x_i, x_step);
658 copy_atm(ctl, atm_i, atm_apr, 0);
659 copy_obs(ctl, obs_i, obs_meas, 0);
660 x2atm(ctl, x_i, atm_i);
661
662 /* Check atmospheric state... */
663 for (int ip = 0; ip < atm_i->np; ip++) {
664 atm_i->p[ip] = GSL_MIN(GSL_MAX(atm_i->p[ip], 5e-7), 5e4);
665 atm_i->t[ip] = GSL_MIN(GSL_MAX(atm_i->t[ip], 100), 400);
666 for (int ig = 0; ig < ctl->ng; ig++)
667 atm_i->q[ig][ip] = GSL_MIN(GSL_MAX(atm_i->q[ig][ip], 0), 1);
668 for (int iw = 0; iw < ctl->nw; iw++)
669 atm_i->k[iw][ip] = GSL_MAX(atm_i->k[iw][ip], 0);
670 }
671
672 /* Forward calculation... */
673 formod(ctl, tbl, atm_i, obs_i);
674 obs2y(ctl, obs_i, y_i, NULL, NULL);
675
676 /* Determine dx = x_i - x_a and dy = y - F(x_i) ... */
677 gsl_vector_memcpy(dx, x_i);
678 gsl_vector_sub(dx, x_a);
679 gsl_vector_memcpy(dy, y_m);
680 gsl_vector_sub(dy, y_i);
681
682 /* Compute cost function... */
683 *chisq = cost_function(dx, dy, s_a_inv, sig_eps_inv);
684
685 /* Modify Levenberg-Marquardt parameter... */
686 if (*chisq > chisq_old) {
687 lmpar *= 10;
688 gsl_vector_sub(x_i, x_step);
689 } else {
690 lmpar /= 10;
691 break;
692 }
693 }
694
695 /* Get normalized step size in state space... */
696 gsl_blas_ddot(x_step, b, &disq);
697 disq /= (double) n;
698
699 /* Convergence test... */
700 if ((it == 1 || it % ret->kernel_recomp == 0) && disq < ret->conv_dmin)
701 break;
702 }
703
704 /* ------------------------------------------------------------
705 Finalize...
706 ------------------------------------------------------------ */
707
708 gsl_matrix_free(a);
709 gsl_matrix_free(cov);
710 gsl_matrix_free(k_i);
711 gsl_matrix_free(s_a_inv);
712
713 gsl_vector_free(b);
714 gsl_vector_free(dx);
715 gsl_vector_free(dy);
716 gsl_vector_free(sig_eps_inv);
717 gsl_vector_free(sig_formod);
718 gsl_vector_free(sig_noise);
719 gsl_vector_free(x_a);
720 gsl_vector_free(x_i);
721 gsl_vector_free(x_step);
722 gsl_vector_free(y_aux);
723 gsl_vector_free(y_i);
724 gsl_vector_free(y_m);
725}
◆ read_nc()
| void read_nc | ( | char * | filename, |
| ncd_t * | ncd | ||
| ) |
Read netCDF file.
Definition at line 729 of file nlte.c.
731 {
732
733 int varid;
734
735 /* Open netCDF file... */
736 printf("Read netCDF file: %s\n", filename);
738
739 /* Read Level-1 data... */
756
757 /* Read Level-2 data... */
764}
double l1_time[L1_NTRACK][L1_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: diff_apr.c:77
◆ main()
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 170 of file nlte.c.
172 {
173
174 static ctl_t ctl;
175 static atm_t atm_apr, atm_clim, atm_i;
176 static obs_t obs_i, obs_meas;
178 static ret_t ret;
179
180 static FILE *in, *out;
181
182 static char filename[LEN], filename2[2 * LEN];
183
185 z[NP];
186
187 static int channel[ND], n, ntask = -1, rank, size;
188
189 /* ------------------------------------------------------------
190 Init...
191 ------------------------------------------------------------ */
192
193 /* MPI... */
194 MPI_Init(&argc, &argv);
195 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
196 MPI_Comm_size(MPI_COMM_WORLD, &size);
197
198 /* Measure CPU time... */
199 TIMER("total", 1);
200
201 /* Check arguments... */
202 if (argc < 3)
203 ERRMSG("Give parameters: <ctl> <filelist>");
204
205 /* Read control parameters... */
206 read_ctl(argc, argv, &ctl);
207 read_ret(argc, argv, &ctl, &ret);
208
209 /* Initialize look-up tables... */
210 tbl_t *tbl = read_tbl(&ctl);
211
212 /* Read retrieval grid... */
213 const int nz = (int) scan_ctl(argc, argv, "NZ", -1, "", NULL);
214 if (nz > NP)
215 ERRMSG("Too many altitudes!");
216 for (int iz = 0; iz < nz; iz++)
217 z[iz] = scan_ctl(argc, argv, "Z", iz, "", NULL);
218
219 /* Read track range... */
220 const int track0 = (int) scan_ctl(argc, argv, "TRACK_MIN", -1, "0", NULL);
221 const int track1 = (int) scan_ctl(argc, argv, "TRACK_MAX", -1, "134", NULL);
222
223 /* Read xtrack range... */
224 const int xtrack0 = (int) scan_ctl(argc, argv, "XTRACK_MIN", -1, "0", NULL);
225 const int xtrack1 =
226 (int) scan_ctl(argc, argv, "XTRACK_MAX", -1, "89", NULL);
227
228 /* Background smoothing... */
229 const double sx = scan_ctl(argc, argv, "SX", -1, "8", NULL);
230 const double sy = scan_ctl(argc, argv, "SY", -1, "2", NULL);
231
232 /* ------------------------------------------------------------
233 Distribute granules...
234 ------------------------------------------------------------ */
235
236 /* Open filelist... */
237 printf("Read filelist: %s\n", argv[2]);
238 if (!(in = fopen(argv[2], "r")))
239 ERRMSG("Cannot open filelist!");
240
241 /* Loop over netCDF files... */
242 while (fscanf(in, "%s", filename) != EOF) {
243
244 /* Distribute files with MPI... */
245 if ((++ntask) % size != rank)
246 continue;
247
248 /* Write info... */
249 printf("Retrieve file %s on rank %d of %d (with %d threads)...\n",
250 filename, rank + 1, size, omp_get_max_threads());
251
252 /* ------------------------------------------------------------
253 Initialize retrieval...
254 ------------------------------------------------------------ */
255
256 /* Read netCDF file... */
257 read_nc(filename, &ncd);
258
259 /* Identify radiance channels... */
260 for (int id = 0; id < ctl.nd; id++) {
261 channel[id] = -999;
264 channel[id] = i;
265 if (channel[id] < 0)
266 ERRMSG("Cannot identify radiance channel!");
267 }
268
269 /* Fill data gaps... */
272
273 /* Set climatological data for center of granule... */
274 atm_clim.np = nz;
275 for (int iz = 0; iz < nz; iz++)
276 atm_clim.z[iz] = z[iz];
277 climatology(&ctl, &atm_clim);
278
279 /* ------------------------------------------------------------
280 Retrieval...
281 ------------------------------------------------------------ */
282
283 /* Loop over swaths... */
284 for (int track = track0; track <= track1; track++) {
285
286 /* Measure CPU time... */
287 TIMER("retrieval", 1);
288
289 /* Loop over scan... */
290 for (int xtrack = xtrack0; xtrack <= xtrack1; xtrack++) {
291
292 /* Store observation data... */
293 obs_meas.nr = 1;
294 obs_meas.time[0] = ncd.l1_time[track][xtrack];
295 obs_meas.obsz[0] = ncd.l1_sat_z[track];
296 obs_meas.obslon[0] = ncd.l1_sat_lon[track];
297 obs_meas.obslat[0] = ncd.l1_sat_lat[track];
298 obs_meas.vplon[0] = ncd.l1_lon[track][xtrack];
299 obs_meas.vplat[0] = ncd.l1_lat[track][xtrack];
300 for (int id = 0; id < ctl.nd; id++)
302
303 /* Flag out 4 micron channels... */
304 for (int id = 0; id < ctl.nd; id++)
305 if (ctl.nu[id] >= 2000)
306 obs_meas.rad[id][0] = GSL_NAN;
307
308 /* Prepare atmospheric data... */
309 copy_atm(&ctl, &atm_apr, &atm_clim, 0);
310 for (int ip = 0; ip < atm_apr.np; ip++) {
311 atm_apr.time[ip] = obs_meas.time[0];
312 atm_apr.lon[ip] = obs_meas.vplon[0];
313 atm_apr.lat[ip] = obs_meas.vplat[0];
314 }
315
316 /* Merge Level-2 data... */
317 init_l2(&ncd, track, xtrack, &ctl, &atm_apr);
318
319 /* Retrieval... */
320 optimal_estimation(&ret, &ctl, tbl, &obs_meas, &obs_i,
321 &atm_apr, &atm_i, &chisq[track][xtrack]);
322
323 /* Run forward model including 4 micron channels... */
324 for (int id = 0; id < ctl.nd; id++)
325 obs_meas.rad[id][0] = obs_i.rad[id][0]
327 formod(&ctl, tbl, &atm_i, &obs_i);
328
329 /* Calculate non-LTE index... */
330 n = 0;
331 ni[track][xtrack] = 0;
332 for (int id = 0; id < ctl.nd; id++)
333 if (ctl.nu[id] >= 2000 && gsl_finite(obs_meas.rad[id][0])) {
334 ni[track][xtrack] +=
337 n++;
338 }
339 ni[track][xtrack] /= n;
340 }
341
342 /* Measure CPU time... */
343 TIMER("retrieval", 3);
344 }
345
346 /* ------------------------------------------------------------
347 Write output...
348 ------------------------------------------------------------ */
349
350 /* Set filename... */
351 sprintf(filename2, "%s.nlte.tab", filename);
352
353 /* Create output file... */
354 printf("Write non-LTE data: %s\n", filename2);
355 if (!(out = fopen(filename2, "w")))
356 ERRMSG("Cannot create file!");
357
358 /* Write header... */
359 fprintf(out,
360 "# $1 = time (seconds since 2000-01-01, 00:00 UTC)\n"
361 "# $2 = longitude [deg]\n"
362 "# $3 = latitude [deg]\n"
363 "# $4 = solar zenith angle [deg]\n"
364 "# $5 = non-LTE index [K]\n" "# $6 = chi^2 of retrieval fit\n");
365
366 /* Write data... */
367 for (int track = track0; track <= track1; track++) {
368 fprintf(out, "\n");
369 for (int xtrack = xtrack0; xtrack <= xtrack1; xtrack++)
370 fprintf(out, "%.2f %g %g %g %g %g\n",
371 ncd.l1_time[track][xtrack],
372 ncd.l1_lon[track][xtrack],
373 ncd.l1_lat[track][xtrack],
375 ncd.l1_lat[track][xtrack]), ni[track][xtrack],
376 chisq[track][xtrack]);
377 }
378
379 /* Close output file... */
380 fclose(out);
381
382 /* Write info... */
383 printf("Retrieval finished on rank %d of %d!\n", rank, size);
384 }
385
386 /* Close file list... */
387 fclose(in);
388
389 /* Free... */
390 free(tbl);
391
392 /* Measure CPU time... */
393 TIMER("total", 3);
394
395 /* Report memory usage... */
396 printf("MEMORY_ATM = %g MByte\n", 4. * sizeof(atm_t) / 1024. / 1024.);
397 printf("MEMORY_CTL = %g MByte\n", 1. * sizeof(ctl_t) / 1024. / 1024.);
399 printf("MEMORY_OBS = %g MByte\n", 3. * sizeof(atm_t) / 1024. / 1024.);
400 printf("MEMORY_RET = %g MByte\n", 1. * sizeof(ret_t) / 1024. / 1024.);
401 printf("MEMORY_TBL = %g MByte\n", 1. * sizeof(tbl_t) / 1024. / 1024.);
402
403 /* Report problem size... */
404 printf("SIZE_TASKS = %d\n", size);
405 printf("SIZE_THREADS = %d\n", omp_get_max_threads());
406
407 /* MPI... */
408 MPI_Finalize();
409
410 return EXIT_SUCCESS;
411}
void fill_gaps(double x[L2_NTRACK][L2_NXTRACK][L2_NLAY], double cx, double cy)
Fill data gaps in L2 data.
Definition: nlte.c:442
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.
Definition: nlte.c:544
void init_l2(ncd_t *ncd, int track, int xtrack, ctl_t *ctl, atm_t *atm)
Initialize with AIRS Level-2 data.
Definition: nlte.c:487
Here is the call graph for this function:
