CrIS Code Collection
Macros | Functions
variance.c File Reference
#include "libcris.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[])
 

Macro Definition Documentation

◆ NLAT_GW

#define NLAT_GW   19

Definition at line 8 of file variance.c.

◆ NLAT_SURF

#define NLAT_SURF   6

Definition at line 9 of file variance.c.

◆ NLAT_TROP

#define NLAT_TROP   73

Definition at line 10 of file variance.c.

◆ NMON

#define NMON   12

Definition at line 13 of file variance.c.

◆ NX

#define NX   3600

Definition at line 16 of file variance.c.

◆ NY

#define NY   1800

Definition at line 19 of file variance.c.

Function Documentation

◆ main()

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

Definition at line 261 of file variance.c.

263 {
264
265 static pert_t *pert;
266
267 static char pertname[LEN], set[LEN];
268
269 const double dc_hlat = 25, dc_tlim = 250;
270
271 static double bt[NX][NY], bt_8mu[NX][NY], bt_8mu_min[NX][NY],
272 bt_8mu_max[NX][NY], dt[NX][NY], mtime[NX][NY], glat[NY], glon[NX],
273 fdc[NX][NY], fwg[NX][NY], fgw[NX][NY], fcw[NX][NY], mean[NX][NY],
274 min[NX][NY], max[NX][NY], var[NX][NY], t_dc, t_gw, help[NX * NY];
275
276 const int iradius = 30, nmin = 10;
277
278 static int n[NX][NY], ndc[NX][NY], ngw[NX][NY], ncw[NX][NY], nwg[NX][NY],
279 det_gw, det_cw, det_dc, det_wg, ilat, ncid, varid, minid, maxid, lonid,
280 latid, npid, dimid[10], help2[NX * NY];
281
282 /* Check arguments... */
283 if (argc < 4)
284 ERRMSG("Give parameters: <ctl> <var.tab> <pert1.nc> [<pert2.nc> ...]");
285
286 /* Get control parameters... */
287 scan_ctl(argc, argv, "SET", -1, "full", set);
288 scan_ctl(argc, argv, "PERTNAME", -1, "4mu", pertname);
289 const int nx = (int) scan_ctl(argc, argv, "NX", -1, "90", NULL);
290 const double lon0 = scan_ctl(argc, argv, "LON0", -1, "-180", NULL);
291 const double lon1 = scan_ctl(argc, argv, "LON1", -1, "180", NULL);
292 const int ny = (int) scan_ctl(argc, argv, "NY", -1, "90", NULL);
293 const double lat0 = scan_ctl(argc, argv, "LAT0", -1, "-90", NULL);
294 const double lat1 = scan_ctl(argc, argv, "LAT1", -1, "90", NULL);
295 const double thresh_gw =
296 scan_ctl(argc, argv, "THRESH_GW", -1, "-999", NULL);
297 const double thresh_dc =
298 scan_ctl(argc, argv, "THRESH_DC", -1, "-999", NULL);
299 const double dt_trop = scan_ctl(argc, argv, "DT_TROP", -1, "0", NULL);
300 const double dt230 = scan_ctl(argc, argv, "DT230", -1, "-999", NULL);
301 const double nu = scan_ctl(argc, argv, "NU", -1, "-999", NULL);
302 const int dc = (int) scan_ctl(argc, argv, "DC", -1, "0", NULL);
303 const int output = (int) scan_ctl(argc, argv, "OUTPUT", -1, "1", NULL);
304
305 /* Allocate... */
306 ALLOC(pert, pert_t, 1);
307
308 /* Check grid dimensions... */
309 if (nx < 1 || nx > NX)
310 ERRMSG("Set 1 <= NX <= MAX!");
311 if (ny < 1 || ny > NY)
312 ERRMSG("Set 1 <= NY <= MAX!");
313
314 /* Loop over perturbation files... */
315 for (int iarg = 3; iarg < argc; iarg++) {
316
317 /* Read perturbation data... */
318 FILE *in;
319 if (!(in = fopen(argv[iarg], "r")))
320 continue;
321 else {
322 fclose(in);
323 read_pert(argv[iarg], pertname, dc, pert);
324 }
325
326 /* Detection... */
327 for (int itrack = 0; itrack < pert->ntrack; itrack++)
328 for (int ixtrack = 0; ixtrack < pert->nxtrack; ixtrack++)
329 for (int ifov = 0; ifov < pert->nfov; ifov++) {
330
331 /* Check data... */
332 if (pert->time[itrack][ixtrack][ifov] < 0
333 || pert->lon[itrack][ixtrack][ifov] < -180
334 || pert->lon[itrack][ixtrack][ifov] > 180
335 || pert->lat[itrack][ixtrack][ifov] < -90
336 || pert->lat[itrack][ixtrack][ifov] > 90
337 || pert->pt[itrack][ixtrack][ifov] < -100
338 || pert->pt[itrack][ixtrack][ifov] > 100
339 || !gsl_finite(pert->bt[itrack][ixtrack][ifov])
340 || !gsl_finite(pert->pt[itrack][ixtrack][ifov])
341 || !gsl_finite(pert->var[itrack][ixtrack][ifov])
342 || !gsl_finite(pert->dc[itrack][ixtrack][ifov]))
343 continue;
344
345 /* Get and check ascending/descending flag... */
346 const int asc =
347 (pert->lat[itrack > 0 ? itrack : itrack + 1][pert->nxtrack / 2]
348 > pert->lat[itrack >
349 0 ? itrack - 1 : itrack][pert->nxtrack / 2]);
350 if (((set[0] == 'a' || set[0] == 'A') && !asc)
351 || ((set[0] == 'd' || set[0] == 'D') && asc))
352 continue;
353
354 /* Check am/pm flag... */
355 const double lt =
356 fmod(pert->time[itrack][ixtrack][ifov], 86400.) / 3600.;
357 if (((set[0] == 'm' || set[0] == 'M') && lt > 12.)
358 || ((set[0] == 'n' || set[0] == 'N') && lt < 12.))
359 continue;
360
361 /* Get grid indices... */
362 int ix = (int) ((pert->lon[itrack][ixtrack][ifov] - lon0)
363 / (lon1 - lon0) * (double) nx);
364 int iy = (int) ((pert->lat[itrack][ixtrack][ifov] - lat0)
365 / (lat1 - lat0) * (double) ny);
366 if (ix < 0 || ix >= nx || iy < 0 || iy >= ny)
367 continue;
368
369 /* Get month index... */
370 int imon =
371 (int) (fmod(pert->time[0][0][0] / 60. / 60. / 24. / 365.25, 1.) *
372 NMON);
373 if (imon < 0 || imon >= NMON)
374 continue;
375
376 /* Get gravity wave detection threshold... */
377 if (thresh_gw <= 0.0) {
378 ilat =
379 locate_irr(t_gw_lat, NLAT_GW, pert->lat[itrack][ixtrack][ifov]);
380 if (asc)
381 t_gw = LIN(t_gw_lat[ilat], t_gw_asc[imon][ilat],
382 t_gw_lat[ilat + 1], t_gw_asc[imon][ilat + 1],
383 pert->lat[itrack][ixtrack][ifov]);
384 else
385 t_gw = LIN(t_gw_lat[ilat], t_gw_dsc[imon][ilat],
386 t_gw_lat[ilat + 1], t_gw_dsc[imon][ilat + 1],
387 pert->lat[itrack][ixtrack][ifov]);
388 } else
389 t_gw = thresh_gw;
390
391 /* Get deep convection detection threshold... */
392 if (thresh_dc <= 0.0) {
393 ilat =
394 locate_irr(t_trop_lat, NLAT_TROP,
395 pert->lat[itrack][ixtrack][ifov]);
396 t_dc =
397 LIN(t_trop_lat[ilat], t_trop[imon][ilat], t_trop_lat[ilat + 1],
398 t_trop[imon][ilat + 1],
399 pert->lat[itrack][ixtrack][ifov]) + dt_trop;
400 } else
401 t_dc = thresh_dc + dt_trop;
402
403 /* Detection of gravity waves... */
404 det_gw = (pert->var[itrack][ixtrack][ifov] >= t_gw);
405
406 /* Detection of convective waves... */
407 det_cw = 0;
408 if (det_gw)
409 for (int itrack2 = GSL_MAX(itrack - iradius, 0);
410 itrack2 <= GSL_MIN(itrack + iradius, pert->ntrack - 1);
411 itrack2++)
412 for (int ixtrack2 = GSL_MAX(ixtrack - iradius, 0);
413 ixtrack2 <= GSL_MIN(ixtrack + iradius, pert->nxtrack - 1);
414 ixtrack2++) {
415 if (det_cw)
416 break;
417 det_cw = (pert->dc[itrack2][ixtrack2][ifov] <= t_dc); // TODO: use ifov2 ?
418 }
419
420 /* Detection of deep convection... */
421 det_dc = (pert->dc[itrack][ixtrack][ifov] <= t_dc);
422
423 /* Detection of wave generation... */
424 det_wg = 0;
425 if (det_dc)
426 for (int itrack2 = GSL_MAX(itrack - iradius, 0);
427 itrack2 <= GSL_MIN(itrack + iradius, pert->ntrack - 1);
428 itrack2++)
429 for (int ixtrack2 = GSL_MAX(ixtrack - iradius, 0);
430 ixtrack2 <= GSL_MIN(ixtrack + iradius, pert->nxtrack - 1);
431 ixtrack2++) {
432 if (det_wg)
433 break;
434 det_wg = (pert->var[itrack2][ixtrack2][ifov] >= t_gw);
435 }
436
437 /* Count events... */
438 n[ix][iy]++;
439 if (det_dc)
440 ndc[ix][iy]++;
441 if (det_wg)
442 nwg[ix][iy]++;
443 if (det_gw)
444 ngw[ix][iy]++;
445 if (det_cw)
446 ncw[ix][iy]++;
447
448 /* Get statistics of perturbations... */
449 mean[ix][iy] += pert->pt[itrack][ixtrack][ifov];
450 var[ix][iy] += gsl_pow_2(pert->pt[itrack][ixtrack][ifov]);
451 max[ix][iy] = GSL_MAX(max[ix][iy], pert->pt[itrack][ixtrack][ifov]);
452 min[ix][iy] = GSL_MIN(min[ix][iy], pert->pt[itrack][ixtrack][ifov]);
453
454 /* Get statistics of brightness temperatures... */
455 bt[ix][iy] += pert->bt[itrack][ixtrack][ifov];
456 bt_8mu[ix][iy] += pert->dc[itrack][ixtrack][ifov];
457 if (n[ix][iy] > 1) {
458 bt_8mu_min[ix][iy] =
459 GSL_MIN(bt_8mu_min[ix][iy], pert->dc[itrack][ixtrack][ifov]);
460 bt_8mu_max[ix][iy] =
461 GSL_MAX(bt_8mu_max[ix][iy], pert->dc[itrack][ixtrack][ifov]);
462 } else {
463 bt_8mu_min[ix][iy] = pert->dc[itrack][ixtrack][ifov];
464 bt_8mu_max[ix][iy] = pert->dc[itrack][ixtrack][ifov];
465 }
466
467 /* Get mean time... */
468 mtime[ix][iy] += pert->time[itrack][ixtrack][ifov];
469 }
470 }
471
472 /* Analyze results... */
473 for (int ix = 0; ix < nx; ix++)
474 for (int iy = 0; iy < ny; iy++) {
475
476 /* Get geolocation... */
477 mtime[ix][iy] /= (double) n[ix][iy];
478 glon[ix]
479 = lon0 + (ix + 0.5) / (double) nx *(
480 lon1 - lon0);
481 glat[iy]
482 = lat0 + (iy + 0.5) / (double) ny *(
483 lat1 - lat0);
484
485 /* Normalize brightness temperatures... */
486 bt[ix][iy] /= (double) n[ix][iy];
487 bt_8mu[ix][iy] /= (double) n[ix][iy];
488
489 /* Get fractions... */
490 fdc[ix][iy] = (double) ndc[ix][iy] / (double) n[ix][iy] * 100.;
491 fwg[ix][iy] = (double) nwg[ix][iy] / (double) ndc[ix][iy] * 100.;
492 fgw[ix][iy] = (double) ngw[ix][iy] / (double) n[ix][iy] * 100.;
493 fcw[ix][iy] = (double) ncw[ix][iy] / (double) ngw[ix][iy] * 100.;
494
495 /* Check number of observations... */
496 if (n[ix][iy] < nmin) {
497 fdc[ix][iy] = GSL_NAN;
498 fwg[ix][iy] = GSL_NAN;
499 fgw[ix][iy] = GSL_NAN;
500 fcw[ix][iy] = GSL_NAN;
501 bt_8mu[ix][iy] = GSL_NAN;
502 bt_8mu_min[ix][iy] = GSL_NAN;
503 bt_8mu_max[ix][iy] = GSL_NAN;
504 }
505
506 /* Check detections of deep convection at high latitudes... */
507 if (fabs(glat[iy]) > dc_hlat && bt_8mu[ix][iy] <= dc_tlim) {
508 fdc[ix][iy] = GSL_NAN;
509 fwg[ix][iy] = GSL_NAN;
510 fcw[ix][iy] = GSL_NAN;
511 }
512
513 /* Estimate noise... */
514 if (dt230 > 0 && nu > 0) {
515 const double nesr = PLANCK(230.0 + dt230, nu) - PLANCK(230.0, nu);
516 dt[ix][iy] = BRIGHT(PLANCK(bt[ix][iy], nu) + nesr, nu) - bt[ix][iy];
517 }
518
519 /* Get mean perturbation and variance... */
520 mean[ix][iy] /= (double) n[ix][iy];
521 var[ix][iy] =
522 var[ix][iy] / (double) n[ix][iy] - gsl_pow_2(mean[ix][iy]);
523 }
524
525 /* Write ASCII file... */
526 if (output == 1) {
527
528 /* Create file... */
529 LOG(1, "Write variance statistics: %s", argv[2]);
530 FILE *out;
531 if (!(out = fopen(argv[2], "w")))
532 ERRMSG("Cannot create file!");
533
534 /* Write header... */
535 fprintf(out,
536 "# $1 = time [s]\n"
537 "# $2 = longitude [deg]\n"
538 "# $3 = latitude [deg]\n"
539 "# $4 = number of footprints\n"
540 "# $5 = fraction of convection events [%%]\n"
541 "# $6 = fraction of wave generating events [%%]\n"
542 "# $7 = fraction of gravity wave events [%%]\n"
543 "# $8 = fraction of convective wave events [%%]\n"
544 "# $9 = mean perturbation [K]\n"
545 "# $10 = minimum perturbation [K]\n");
546 fprintf(out,
547 "# $11 = maximum perturbation [K]\n"
548 "# $12 = variance [K^2]\n"
549 "# $13 = mean surface temperature [K]\n"
550 "# $14 = minimum surface temperature [K]\n"
551 "# $15 = maximum surface temperature [K]\n"
552 "# $16 = mean background temperature [K]\n"
553 "# $17 = noise estimate [K]\n");
554
555 /* Write results... */
556 for (int iy = 0; iy < ny; iy++) {
557 if (iy == 0 || nx > 1)
558 fprintf(out, "\n");
559 for (int ix = 0; ix < nx; ix++)
560 fprintf(out, "%.2f %g %g %d %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
561 mtime[ix][iy], glon[ix], glat[iy], n[ix][iy],
562 fdc[ix][iy], fwg[ix][iy], fgw[ix][iy], fcw[ix][iy],
563 mean[ix][iy], min[ix][iy], max[ix][iy], var[ix][iy],
564 bt_8mu[ix][iy], bt_8mu_min[ix][iy], bt_8mu_max[ix][iy],
565 bt[ix][iy], dt[ix][iy]);
566 }
567
568 /* Close file... */
569 fclose(out);
570 }
571
572 /* Write netCDF file... */
573 else if (output == 2) {
574
575 /* Create netCDF file... */
576 LOG(1, "Write variance statistics: %s", argv[2]);
577 NC(nc_create(argv[2], NC_CLOBBER, &ncid));
578
579 /* Set dimensions... */
580 NC(nc_def_dim(ncid, "lat", (size_t) ny, &dimid[0]));
581 NC(nc_def_dim(ncid, "lon", (size_t) nx, &dimid[1]));
582
583 /* Add variables... */
584 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 1, &dimid[0], &latid));
585 add_att(ncid, latid, "deg", "latitude");
586 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 1, &dimid[1], &lonid));
587 add_att(ncid, lonid, "deg", "longitude");
588 NC(nc_def_var(ncid, "var", NC_FLOAT, 2, dimid, &varid));
589 add_att(ncid, varid, "K^2", "brightness temperature variance");
590 NC(nc_def_var(ncid, "min", NC_FLOAT, 2, dimid, &minid));
591 add_att(ncid, minid, "K", "brightness temperature minimum");
592 NC(nc_def_var(ncid, "max", NC_FLOAT, 2, dimid, &maxid));
593 add_att(ncid, maxid, "K", "brightness temperature maximum");
594 NC(nc_def_var(ncid, "np", NC_INT, 2, dimid, &npid));
595 add_att(ncid, npid, "1", "number of footprints");
596
597 /* Leave define mode... */
598 NC(nc_enddef(ncid));
599
600 /* Write data... */
601 NC(nc_put_var_double(ncid, latid, glat));
602 NC(nc_put_var_double(ncid, lonid, glon));
603 for (int ix = 0; ix < nx; ix++)
604 for (int iy = 0; iy < ny; iy++)
605 help[iy * nx + ix] = var[ix][iy] - POW2(dt[ix][iy]);
606 NC(nc_put_var_double(ncid, varid, help));
607 for (int ix = 0; ix < nx; ix++)
608 for (int iy = 0; iy < ny; iy++)
609 help[iy * nx + ix] = min[ix][iy];
610 NC(nc_put_var_double(ncid, minid, help));
611 for (int ix = 0; ix < nx; ix++)
612 for (int iy = 0; iy < ny; iy++)
613 help[iy * nx + ix] = max[ix][iy];
614 NC(nc_put_var_double(ncid, maxid, help));
615 for (int ix = 0; ix < nx; ix++)
616 for (int iy = 0; iy < ny; iy++)
617 help2[iy * nx + ix] = n[ix][iy];
618 NC(nc_put_var_int(ncid, npid, help2));
619
620 /* Close file... */
621 NC(nc_close(ncid));
622 }
623
624 else
625 ERRMSG("Unknown output format!");
626
627 /* Free... */
628 free(pert);
629
630 return EXIT_SUCCESS;
631}
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
LOG
#define LOG(level,...)
Print log message.
Definition: jurassic.h:221
LIN
#define LIN(x0, y0, x1, y1, x)
Compute linear interpolation.
Definition: jurassic.h:164
add_att
void add_att(const int ncid, const int varid, const char *unit, const char *long_name)
Add variable attributes to netCDF file.
Definition: libcris.c:5
read_pert
void read_pert(char *filename, char *pertname, int dc, pert_t *pert)
Read radiance perturbation data.
Definition: libcris.c:1270
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libcris.h:61
pert_t
Perturbation data.
Definition: libcris.h:131
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature (cloud channel) [K].
Definition: libcris.h:152
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Latitude [deg].
Definition: libcris.h:149
pert_t::nfov
int nfov
Number of field of views.
Definition: libcris.h:140
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libcris.h:158
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libcris.h:134
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libcris.h:161
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Time (seconds since 2000-01-01T00:00Z).
Definition: libcris.h:143
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature (4.3 or 15 micron) [K].
Definition: libcris.h:155
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libcris.h:137
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Longitude [deg].
Definition: libcris.h:146
NLAT_GW
#define NLAT_GW
Definition: variance.c:8
NLAT_TROP
#define NLAT_TROP
Definition: variance.c:10
NX
#define NX
Definition: variance.c:16
NY
#define NY
Definition: variance.c:19
NMON
#define NMON
Definition: variance.c:13
Here is the call graph for this function:
Generated on Wed Feb 5 2025 20:00:07 for CrIS Code Collection by doxygen 1.9.4