AIRS Code Collection
Macros | Functions
variance.c File Reference

Calculate brightness temperature variances. More...

#include "libairs.h"

Go to the source code of this file.

Macros

#define NLAT_GW   19
 
#define NLAT_SURF   6
 
#define NLAT_TROP   73
 
#define NMON   12
 
#define NX   3600
 
#define NY   1800
 

Functions

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

Detailed Description

Calculate brightness temperature variances.

Definition in file variance.c.

Macro Definition Documentation

◆ NLAT_GW

#define NLAT_GW   19

Definition at line 33 of file variance.c.

◆ NLAT_SURF

#define NLAT_SURF   6

Definition at line 34 of file variance.c.

◆ NLAT_TROP

#define NLAT_TROP   73

Definition at line 35 of file variance.c.

◆ NMON

#define NMON   12

Definition at line 38 of file variance.c.

◆ NX

#define NX   3600

Definition at line 41 of file variance.c.

◆ NY

#define NY   1800

Definition at line 44 of file variance.c.

Function Documentation

◆ main()

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

Definition at line 286 of file variance.c.

288 {
289
290 static pert_t *pert;
291
292 static wave_t *wave;
293
294 static char pertname[LEN], set[LEN];
295
296 const double dc_hlat = 25, dc_tlim = 250;
297
298 static double bt[NX][NY], bt_8mu[NX][NY], bt_8mu_min[NX][NY],
299 bt_8mu_max[NX][NY], dt[NX][NY], mtime[NX][NY], glat[NY], glon[NX],
300 fdc[NX][NY], fwg[NX][NY], fgw[NX][NY], fcw[NX][NY], mean[NX][NY],
301 min[NX][NY], max[NX][NY], var[NX][NY], t_dc, t_gw, help[NX * NY];
302
303 const int nmin = 10, iradius = 30;
304
305 static int n[NX][NY], ndc[NX][NY], ngw[NX][NY], ncw[NX][NY], nwg[NX][NY],
306 det_gw, det_cw, det_dc, det_wg, ilat, ncid, varid, minid, maxid, lonid,
307 latid, npid, dimid[10], help2[NX * NY];
308
309 /* Check arguments... */
310 if (argc < 4)
311 ERRMSG("Give parameters: <ctl> <var.tab> <pert1.nc> [<pert2.nc> ...]");
312
313 /* Get control parameters... */
314 scan_ctl(argc, argv, "SET", -1, "full", set);
315 scan_ctl(argc, argv, "PERTNAME", -1, "4mu", pertname);
316 const int nx = (int) scan_ctl(argc, argv, "NX", -1, "360", NULL);
317 const double lon0 = scan_ctl(argc, argv, "LON0", -1, "-180", NULL);
318 const double lon1 = scan_ctl(argc, argv, "LON1", -1, "180", NULL);
319 const int ny = (int) scan_ctl(argc, argv, "NY", -1, "180", NULL);
320 const double lat0 = scan_ctl(argc, argv, "LAT0", -1, "-90", NULL);
321 const double lat1 = scan_ctl(argc, argv, "LAT1", -1, "90", NULL);
322 const int bg_poly_x =
323 (int) scan_ctl(argc, argv, "BG_POLY_X", -1, "0", NULL);
324 const int bg_poly_y =
325 (int) scan_ctl(argc, argv, "BG_POLY_Y", -1, "0", NULL);
326 const int bg_smooth_x =
327 (int) scan_ctl(argc, argv, "BG_SMOOTH_X", -1, "0", NULL);
328 const int bg_smooth_y =
329 (int) scan_ctl(argc, argv, "BG_SMOOTH_Y", -1, "0", NULL);
330 const double gauss_fwhm = scan_ctl(argc, argv, "GAUSS_FWHM", -1, "0", NULL);
331 const double var_dh = scan_ctl(argc, argv, "VAR_DH", -1, "0", NULL);
332 const double thresh_gw =
333 scan_ctl(argc, argv, "THRESH_GW", -1, "-999", NULL);
334 const double thresh_dc =
335 scan_ctl(argc, argv, "THRESH_DC", -1, "-999", NULL);
336 const double dt_trop = scan_ctl(argc, argv, "DT_TROP", -1, "0", NULL);
337 const double dt230 = scan_ctl(argc, argv, "DT230", -1, "0.16", NULL);
338 const double nu = scan_ctl(argc, argv, "NU", -1, "2345.0", NULL);
339 const int output = (int) scan_ctl(argc, argv, "OUTPUT", -1, "1", NULL);
340
341 /* Allocate... */
342 ALLOC(pert, pert_t, 1);
343
344 /* Check grid dimensions... */
345 if (nx < 1 || nx > NX)
346 ERRMSG("Set 1 <= NX <= MAX!");
347 if (ny < 1 || ny > NY)
348 ERRMSG("Set 1 <= NY <= MAX!");
349
350 /* Loop over perturbation files... */
351 for (int iarg = 3; iarg < argc; iarg++) {
352
353 /* Read perturbation data... */
354 FILE *in;
355 if (!(in = fopen(argv[iarg], "r")))
356 continue;
357 else {
358 fclose(in);
359 read_pert(argv[iarg], pertname, pert);
360 }
361
362 /* Recalculate background and perturbations... */
363 if (bg_poly_x > 0 || bg_poly_y > 0 ||
364 bg_smooth_x > 0 || bg_smooth_y > 0 || gauss_fwhm > 0 || var_dh > 0) {
365
366 /* Allocate... */
367 ALLOC(wave, wave_t, 1);
368
369 /* Convert to wave analysis struct... */
370 pert2wave(pert, wave, 0, pert->ntrack - 1, 0, pert->nxtrack - 1);
371
372 /* Estimate background... */
373 background_poly(wave, bg_poly_x, bg_poly_y);
374 background_smooth(wave, bg_smooth_x, bg_smooth_y);
375
376 /* Gaussian filter... */
377 gauss(wave, gauss_fwhm);
378
379 /* Compute variance... */
380 variance(wave, var_dh);
381
382 /* Copy data... */
383 for (int ix = 0; ix < wave->nx; ix++)
384 for (int iy = 0; iy < wave->ny; iy++) {
385 pert->pt[iy][ix] = wave->pt[ix][iy];
386 pert->var[iy][ix] = wave->var[ix][iy];
387 }
388
389 /* Free... */
390 free(wave);
391 }
392
393 /* Detection... */
394 for (int itrack = 0; itrack < pert->ntrack; itrack++)
395 for (int ixtrack = 0; ixtrack < pert->nxtrack; ixtrack++) {
396
397 /* Check data... */
398 if (pert->time[itrack][ixtrack] < 0
399 || pert->lon[itrack][ixtrack] < -180
400 || pert->lon[itrack][ixtrack] > 180
401 || pert->lat[itrack][ixtrack] < -90
402 || pert->lat[itrack][ixtrack] > 90
403 || pert->pt[itrack][ixtrack] < -100
404 || pert->pt[itrack][ixtrack] > 100
405 || !gsl_finite(pert->bt[itrack][ixtrack])
406 || !gsl_finite(pert->pt[itrack][ixtrack])
407 || !gsl_finite(pert->var[itrack][ixtrack])
408 || !gsl_finite(pert->dc[itrack][ixtrack]))
409 continue;
410
411 /* Get and check ascending/descending flag... */
412 int asc =
413 (pert->lat[itrack > 0 ? itrack : itrack + 1][pert->nxtrack / 2] >
414 pert->lat[itrack > 0 ? itrack - 1 : itrack][pert->nxtrack / 2]);
415 if (((set[0] == 'a' || set[0] == 'A') && !asc)
416 || ((set[0] == 'd' || set[0] == 'D') && asc))
417 continue;
418
419 /* Check am/pm flag... */
420 double lt = fmod(pert->time[itrack][ixtrack], 86400.) / 3600.;
421 if (((set[0] == 'm' || set[0] == 'M') && lt > 12.)
422 || ((set[0] == 'n' || set[0] == 'N') && lt < 12.))
423 continue;
424
425 /* Get grid indices... */
426 int ix =
427 (int) ((pert->lon[itrack][ixtrack] - lon0) / (lon1 -
428 lon0) * (double) nx);
429 int iy =
430 (int) ((pert->lat[itrack][ixtrack] - lat0) / (lat1 -
431 lat0) * (double) ny);
432 if (ix < 0 || ix >= nx || iy < 0 || iy >= ny)
433 continue;
434
435 /* Get month index... */
436 int imon =
437 (int) (fmod(pert->time[0][0] / 60. / 60. / 24. / 365.25, 1.) *
438 NMON);
439 if (imon < 0 || imon >= NMON)
440 continue;
441
442 /* Get gravity wave detection threshold... */
443 if (thresh_gw <= 0.0) {
444 ilat = locate_irr(t_gw_lat, NLAT_GW, pert->lat[itrack][ixtrack]);
445 if (asc)
446 t_gw = LIN(t_gw_lat[ilat], t_gw_asc[imon][ilat],
447 t_gw_lat[ilat + 1], t_gw_asc[imon][ilat + 1],
448 pert->lat[itrack][ixtrack]);
449 else
450 t_gw = LIN(t_gw_lat[ilat], t_gw_dsc[imon][ilat],
451 t_gw_lat[ilat + 1], t_gw_dsc[imon][ilat + 1],
452 pert->lat[itrack][ixtrack]);
453 } else
454 t_gw = thresh_gw;
455
456 /* Get deep convection detection threshold... */
457 if (thresh_dc <= 0.0) {
458 ilat =
459 locate_irr(t_trop_lat, NLAT_TROP, pert->lat[itrack][ixtrack]);
460 t_dc =
461 LIN(t_trop_lat[ilat], t_trop[imon][ilat], t_trop_lat[ilat + 1],
462 t_trop[imon][ilat + 1], pert->lat[itrack][ixtrack]) + dt_trop;
463 } else
464 t_dc = thresh_dc + dt_trop;
465
466 /* Detection of gravity waves... */
467 det_gw = (pert->var[itrack][ixtrack] >= t_gw);
468
469 /* Detection of convective waves... */
470 det_cw = 0;
471 if (det_gw)
472 for (int itrack2 = GSL_MAX(itrack - iradius, 0);
473 itrack2 <= GSL_MIN(itrack + iradius, pert->ntrack - 1);
474 itrack2++)
475 for (int ixtrack2 = GSL_MAX(ixtrack - iradius, 0);
476 ixtrack2 <= GSL_MIN(ixtrack + iradius, pert->nxtrack - 1);
477 ixtrack2++) {
478 if (det_cw)
479 break;
480 det_cw = (pert->dc[itrack2][ixtrack2] <= t_dc);
481 }
482
483 /* Detection of deep convection... */
484 det_dc = (pert->dc[itrack][ixtrack] <= t_dc);
485
486 /* Detection of wave generation... */
487 det_wg = 0;
488 if (det_dc)
489 for (int itrack2 = GSL_MAX(itrack - iradius, 0);
490 itrack2 <= GSL_MIN(itrack + iradius, pert->ntrack - 1);
491 itrack2++)
492 for (int ixtrack2 = GSL_MAX(ixtrack - iradius, 0);
493 ixtrack2 <= GSL_MIN(ixtrack + iradius, pert->nxtrack - 1);
494 ixtrack2++) {
495 if (det_wg)
496 break;
497 det_wg = (pert->var[itrack2][ixtrack2] >= t_gw);
498 }
499
500 /* Count events... */
501 n[ix][iy]++;
502 if (det_dc)
503 ndc[ix][iy]++;
504 if (det_wg)
505 nwg[ix][iy]++;
506 if (det_gw)
507 ngw[ix][iy]++;
508 if (det_cw)
509 ncw[ix][iy]++;
510
511 /* Get statistics of perturbations... */
512 mean[ix][iy] += pert->pt[itrack][ixtrack];
513 var[ix][iy] += gsl_pow_2(pert->pt[itrack][ixtrack]);
514 max[ix][iy] = GSL_MAX(max[ix][iy], pert->pt[itrack][ixtrack]);
515 min[ix][iy] = GSL_MIN(min[ix][iy], pert->pt[itrack][ixtrack]);
516
517 /* Get statistics of brightness temperatures... */
518 bt[ix][iy] += pert->bt[itrack][ixtrack];
519 bt_8mu[ix][iy] += pert->dc[itrack][ixtrack];
520 if (n[ix][iy] > 1) {
521 bt_8mu_min[ix][iy]
522 = GSL_MIN(bt_8mu_min[ix][iy], pert->dc[itrack][ixtrack]);
523 bt_8mu_max[ix][iy]
524 = GSL_MAX(bt_8mu_max[ix][iy], pert->dc[itrack][ixtrack]);
525 } else {
526 bt_8mu_min[ix][iy] = pert->dc[itrack][ixtrack];
527 bt_8mu_max[ix][iy] = pert->dc[itrack][ixtrack];
528 }
529
530 /* Get mean time... */
531 mtime[ix][iy] += pert->time[itrack][ixtrack];
532 }
533 }
534
535 /* Analyze results... */
536 for (int ix = 0; ix < nx; ix++)
537 for (int iy = 0; iy < ny; iy++) {
538
539 /* Get geolocation... */
540 mtime[ix][iy] /= (double) n[ix][iy];
541 glon[ix]
542 = lon0 + (ix + 0.5) / (double) nx *(
543 lon1 - lon0);
544 glat[iy]
545 = lat0 + (iy + 0.5) / (double) ny *(
546 lat1 - lat0);
547
548 /* Normalize brightness temperatures... */
549 bt[ix][iy] /= (double) n[ix][iy];
550 bt_8mu[ix][iy] /= (double) n[ix][iy];
551
552 /* Get fractions... */
553 fdc[ix][iy] = (double) ndc[ix][iy] / (double) n[ix][iy] * 100.;
554 fwg[ix][iy] = (double) nwg[ix][iy] / (double) ndc[ix][iy] * 100.;
555 fgw[ix][iy] = (double) ngw[ix][iy] / (double) n[ix][iy] * 100.;
556 fcw[ix][iy] = (double) ncw[ix][iy] / (double) ngw[ix][iy] * 100.;
557
558 /* Check number of observations... */
559 if (n[ix][iy] < nmin) {
560 fdc[ix][iy] = GSL_NAN;
561 fwg[ix][iy] = GSL_NAN;
562 fgw[ix][iy] = GSL_NAN;
563 fcw[ix][iy] = GSL_NAN;
564 bt_8mu[ix][iy] = GSL_NAN;
565 bt_8mu_min[ix][iy] = GSL_NAN;
566 bt_8mu_max[ix][iy] = GSL_NAN;
567 }
568
569 /* Check detections of deep convection at high latitudes... */
570 if (fabs(glat[iy]) > dc_hlat && bt_8mu[ix][iy] <= dc_tlim) {
571 fdc[ix][iy] = GSL_NAN;
572 fwg[ix][iy] = GSL_NAN;
573 fcw[ix][iy] = GSL_NAN;
574 }
575
576 /* Estimate noise... */
577 if (dt230 > 0) {
578 const double nesr = PLANCK(230.0 + dt230, nu) - PLANCK(230.0, nu);
579 dt[ix][iy] = BRIGHT(PLANCK(bt[ix][iy], nu) + nesr, nu) - bt[ix][iy];
580 }
581
582 /* Get mean perturbation and variance... */
583 mean[ix][iy] /= (double) n[ix][iy];
584 var[ix][iy] =
585 var[ix][iy] / (double) n[ix][iy] - gsl_pow_2(mean[ix][iy]);
586 }
587
588 /* Write ASCII file... */
589 if (output == 1) {
590
591 /* Write info... */
592 printf("Write variance statistics: %s\n", argv[2]);
593
594 /* Create file... */
595 FILE *out;
596 if (!(out = fopen(argv[2], "w")))
597 ERRMSG("Cannot create file!");
598
599 /* Write header... */
600 fprintf(out,
601 "# $1 = time [s]\n"
602 "# $2 = longitude [deg]\n"
603 "# $3 = latitude [deg]\n"
604 "# $4 = number of footprints\n"
605 "# $5 = fraction of convection events [%%]\n"
606 "# $6 = fraction of wave generating events [%%]\n"
607 "# $7 = fraction of gravity wave events [%%]\n"
608 "# $8 = fraction of convective wave events [%%]\n"
609 "# $9 = mean perturbation [K]\n"
610 "# $10 = minimum perturbation [K]\n");
611 fprintf(out,
612 "# $11 = maximum perturbation [K]\n"
613 "# $12 = variance [K^2]\n"
614 "# $13 = mean surface temperature [K]\n"
615 "# $14 = minimum surface temperature [K]\n"
616 "# $15 = maximum surface temperature [K]\n"
617 "# $16 = mean background temperature [K]\n"
618 "# $17 = noise estimate [K]\n");
619
620 /* Write results... */
621 for (int iy = 0; iy < ny; iy++) {
622 if (iy == 0 || nx > 1)
623 fprintf(out, "\n");
624 for (int ix = 0; ix < nx; ix++)
625 fprintf(out, "%.2f %g %g %d %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
626 mtime[ix][iy], glon[ix], glat[iy], n[ix][iy],
627 fdc[ix][iy], fwg[ix][iy], fgw[ix][iy], fcw[ix][iy],
628 mean[ix][iy], min[ix][iy], max[ix][iy], var[ix][iy],
629 bt_8mu[ix][iy], bt_8mu_min[ix][iy], bt_8mu_max[ix][iy],
630 bt[ix][iy], dt[ix][iy]);
631 }
632
633 /* Close file... */
634 fclose(out);
635 }
636
637 /* Write netCDF file... */
638 else if (output == 2) {
639
640 /* Create netCDF file... */
641 printf("Write variance statistics: %s\n", argv[2]);
642 NC(nc_create(argv[2], NC_CLOBBER, &ncid));
643
644 /* Set dimensions... */
645 NC(nc_def_dim(ncid, "lat", (size_t) ny, &dimid[0]));
646 NC(nc_def_dim(ncid, "lon", (size_t) nx, &dimid[1]));
647
648 /* Add variables... */
649 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 1, &dimid[0], &latid));
650 add_att(ncid, latid, "deg", "latitude");
651 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 1, &dimid[1], &lonid));
652 add_att(ncid, lonid, "deg", "longitude");
653 NC(nc_def_var(ncid, "var", NC_FLOAT, 2, dimid, &varid));
654 add_att(ncid, varid, "K^2", "brightness temperature variance");
655 NC(nc_def_var(ncid, "min", NC_FLOAT, 2, dimid, &minid));
656 add_att(ncid, minid, "K", "brightness temperature minimum");
657 NC(nc_def_var(ncid, "max", NC_FLOAT, 2, dimid, &maxid));
658 add_att(ncid, maxid, "K", "brightness temperature maximum");
659 NC(nc_def_var(ncid, "np", NC_INT, 2, dimid, &npid));
660 add_att(ncid, npid, "1", "number of footprints");
661
662 /* Leave define mode... */
663 NC(nc_enddef(ncid));
664
665 /* Write data... */
666 NC(nc_put_var_double(ncid, latid, glat));
667 NC(nc_put_var_double(ncid, lonid, glon));
668 for (int ix = 0; ix < nx; ix++)
669 for (int iy = 0; iy < ny; iy++)
670 help[iy * nx + ix] = var[ix][iy] - POW2(dt[ix][iy]);
671 NC(nc_put_var_double(ncid, varid, help));
672 for (int ix = 0; ix < nx; ix++)
673 for (int iy = 0; iy < ny; iy++)
674 help[iy * nx + ix] = min[ix][iy];
675 NC(nc_put_var_double(ncid, minid, help));
676 for (int ix = 0; ix < nx; ix++)
677 for (int iy = 0; iy < ny; iy++)
678 help[iy * nx + ix] = max[ix][iy];
679 NC(nc_put_var_double(ncid, maxid, help));
680 for (int ix = 0; ix < nx; ix++)
681 for (int iy = 0; iy < ny; iy++)
682 help2[iy * nx + ix] = n[ix][iy];
683 NC(nc_put_var_int(ncid, npid, help2));
684
685 /* Close file... */
686 NC(nc_close(ncid));
687 }
688
689 else
690 ERRMSG("Unknown output format!");
691
692 /* Free... */
693 free(pert);
694
695 return EXIT_SUCCESS;
696}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: diff_apr.c:35
locate_irr
int locate_irr(const double *xx, const int n, const double x)
Find array index for irregular grid.
Definition: jurassic.c:4103
scan_ctl
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
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:383
BRIGHT
#define BRIGHT(rad, nu)
Compute brightness temperature.
Definition: jurassic.h:126
POW2
#define POW2(x)
Compute x^2.
Definition: jurassic.h:187
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory.
Definition: jurassic.h:121
PLANCK
#define PLANCK(T, nu)
Compute Planck function.
Definition: jurassic.h:183
LIN
#define LIN(x0, y0, x1, y1, x)
Compute linear interpolation.
Definition: jurassic.h:164
background_smooth
void background_smooth(wave_t *wave, int npts_x, int npts_y)
Smooth background.
Definition: libairs.c:176
add_att
void add_att(int ncid, int varid, const char *unit, const char *long_name)
Add variable attributes to netCDF file.
Definition: libairs.c:30
background_poly
void background_poly(wave_t *wave, int dim_x, int dim_y)
Get background based on polynomial fits.
Definition: libairs.c:126
variance
void variance(wave_t *wave, double dh)
Compute local variance.
Definition: libairs.c:1584
pert2wave
void pert2wave(pert_t *pert, wave_t *wave, int track0, int track1, int xtrack0, int xtrack1)
Convert radiance perturbation data to wave analysis struct.
Definition: libairs.c:999
read_pert
void read_pert(char *filename, char *pertname, pert_t *pert)
Read radiance perturbation data.
Definition: libairs.c:1137
gauss
void gauss(wave_t *wave, double fwhm)
Apply Gaussian filter to perturbations...
Definition: libairs.c:579
pert_t
Perturbation data.
Definition: libairs.h:205
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libairs.h:232
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libairs.h:229
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libairs.h:214
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libairs.h:223
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libairs.h:208
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libairs.h:211
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK]
Longitude [deg].
Definition: libairs.h:217
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (4 or 15 micron) [K].
Definition: libairs.h:226
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK]
Latitude [deg].
Definition: libairs.h:220
wave_t
Wave analysis data.
Definition: libairs.h:287
wave_t::nx
int nx
Number of across-track values.
Definition: libairs.h:290
wave_t::ny
int ny
Number of along-track values.
Definition: libairs.h:293
wave_t::var
double var[WX][WY]
Variance [K].
Definition: libairs.h:323
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libairs.h:320
NLAT_GW
#define NLAT_GW
Definition: variance.c:33
NLAT_TROP
#define NLAT_TROP
Definition: variance.c:35
NX
#define NX
Definition: variance.c:41
NY
#define NY
Definition: variance.c:44
NMON
#define NMON
Definition: variance.c:38
Here is the call graph for this function:
Generated on Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4