CrIS Code Collection
Macros | Functions
perturbation.c File Reference
#include "libcris.h"

Go to the source code of this file.

Macros

#define N4   54
 Number of 4.3 micron channels: More...
 
#define N15_LOW   15
 Number of 15 micron low channels: More...
 
#define N15_HIGH   1
 Number of 15 micron high channels: More...
 
#define N8   1
 Number of 8.1 micron channels: More...
 
#define N10   1
 Number of 10.4 micron channels: More...
 
#define LOOP_ALL(ntrack)
 Loop over all footprints. More...
 
#define WRITE_FIELD(expr, varid)
 
#define OUTLIER_FILTER(BT, DTMED)
 
#define APPLY_BIAS_CORR(pert)
 

Functions

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

Macro Definition Documentation

◆ N4

#define N4   54

Number of 4.3 micron channels:

Definition at line 8 of file perturbation.c.

◆ N15_LOW

#define N15_LOW   15

Number of 15 micron low channels:

Definition at line 11 of file perturbation.c.

◆ N15_HIGH

#define N15_HIGH   1

Number of 15 micron high channels:

Definition at line 14 of file perturbation.c.

◆ N8

#define N8   1

Number of 8.1 micron channels:

Definition at line 17 of file perturbation.c.

◆ N10

#define N10   1

Number of 10.4 micron channels:

Definition at line 20 of file perturbation.c.

◆ LOOP_ALL

#define LOOP_ALL (   ntrack)
Value:
for (int track = 0; track < ntrack; track++) \
for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++) \
for (int ifov = 0; ifov < L1_NFOV; ifov++)
L1_NXTRACK
#define L1_NXTRACK
Across-track size of CrIS radiance granule.
Definition: libcris.h:24
L1_NFOV
#define L1_NFOV
Number of field of views of CrIS radiance granule.
Definition: libcris.h:27

Loop over all footprints.

Definition at line 27 of file perturbation.c.

◆ WRITE_FIELD

#define WRITE_FIELD (   expr,
  varid 
)
Value:
do { \
n = 0; \
LOOP_ALL(pert_4mu->ntrack) \
help[n++] = (expr); \
NC(nc_put_var_double(ncid, (varid), help)); \
} while (0)

Definition at line 33 of file perturbation.c.

◆ OUTLIER_FILTER

#define OUTLIER_FILTER (   BT,
  DTMED 
)

◆ APPLY_BIAS_CORR

#define APPLY_BIAS_CORR (   pert)
Value:
do { \
LOOP_ALL((pert)->ntrack) { \
const int idx = (track * L1_NXTRACK + xtrack) * L1_NFOV + ifov; \
help[idx] = 0.0; \
double wsum = 0.0; \
for (int dtrack2 = -radius; dtrack2 <= radius; ++dtrack2) { \
int t2 = track + dtrack2; \
if (t2 < 0 || t2 >= (pert)->ntrack) \
continue; \
if (!isfinite((pert)->pt[t2][xtrack][ifov])) \
continue; \
const double w = exp(-(double)(dtrack2 * dtrack2) * inv2sigma2); \
help[idx] += w * (pert)->pt[t2][xtrack][ifov]; \
wsum += w; \
} \
if (wsum > 0.0) help[idx] /= wsum; \
} \
LOOP_ALL((pert)->ntrack) { \
const int idx = (track * L1_NXTRACK + xtrack) * L1_NFOV + ifov; \
(pert)->pt[track][xtrack][ifov] -= help[idx]; \
} \
} while (0)
LOOP_ALL
#define LOOP_ALL(ntrack)
Loop over all footprints.
Definition: perturbation.c:27

Function Documentation

◆ main()

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

Definition at line 44 of file perturbation.c.

46 {
47
48 static cris_l1_t l1;
49
50 static pert_t *pert_4mu, *pert_15mu_low, *pert_15mu_high;
51
52 const double var_dh2 = 100. * 100.;
53
54 static double help[PERT_NTRACK * PERT_NXTRACK * PERT_NFOV], rad, nu,
55 x[L1_NXTRACK], y[L1_NXTRACK], x0[3], x1[3];
56
57 const int list_4mu[N4]
58 = { 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283,
59 284, 285, 286, 288, 289, 291, 292, 294, 295, 296, 297, 298, 299, 300,
60 302, 303, 304, 305, 306, 307, 318, 319, 321, 322, 323, 324, 325, 326,
61 328, 330, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341
62 };
63
64 const int list_15mu_low[N15_LOW]
65 = { 5, 8, 10, 13, 15, 17, 23, 34, 36, 39, 41, 44, 46, 49, 51 };
66
67 const int list_15mu_high[N15_HIGH] = { 31 };
68
69 const int list_8mu[N8] = { 36 };
70
71 const int list_10mu[N10] = { 502 };
72
73 const int dtrack = 3, dxtrack = 3;
74
75 static int dimid[3], n, ncid, track0,
76 time_varid, lon_varid, lat_varid, bt_4mu_varid, init,
77 bt_4mu_pt_varid, bt_4mu_var_varid, bt_8mu_varid, bt_10mu_varid,
78 bt_15mu_low_varid, bt_15mu_low_pt_varid, bt_15mu_low_var_varid,
79 bt_15mu_high_varid, bt_15mu_high_pt_varid, bt_15mu_high_var_varid;
80
81 /* Check arguments... */
82 if (argc < 4)
83 ERRMSG("Give parameters: <ctl> <pert.nc> <l1b_file1> [<l1b_file2> ...]");
84
85 /* Get control parameters... */
86 const int apo = (int) scan_ctl(argc, argv, "APO", -1, "0", NULL);
87 const int bias = (int) scan_ctl(argc, argv, "BIAS", -1, "1", NULL);
88 const double dtmed = scan_ctl(argc, argv, "DTMED", -1, "10.0", NULL);
89
90 /* Allocate... */
91 ALLOC(pert_4mu, pert_t, 1);
92 ALLOC(pert_15mu_low, pert_t, 1);
93 ALLOC(pert_15mu_high, pert_t, 1);
94
95 /* ------------------------------------------------------------
96 Read CrIS Level-1B files...
97 ------------------------------------------------------------ */
98
99 /* Loop over files... */
100 for (int iarg = 3; iarg < argc; iarg++) {
101
102 /* Read CrIS data... */
103 if (!read_cris_l1(argv[iarg], &l1, apo))
104 continue;
105
106 /* Write channel information... */
107 if (!init) {
108 init = 1;
109
110 /* 4.3 micron channels... */
111 LOG(2, "4.3 micron channels:");
112 n = 0;
113 nu = 0;
114 for (int i = 0; i < N4; i++) {
115 nu += l1.nu_sw[list_4mu[i]];
116 n++;
117 LOG(2, " count= %2d | channel index= SW_%03d | nu= %.4f cm^-1", n,
118 list_4mu[i], l1.nu_sw[list_4mu[i]]);
119 }
120 LOG(2, " nu_mean= %.4f cm^-1", nu / n);
121
122 /* 15 micron low channels... */
123 LOG(2, "15 micron low channels:");
124 n = 0;
125 nu = 0;
126 for (int i = 0; i < N15_LOW; i++) {
127 nu += l1.nu_lw[list_15mu_low[i]];
128 n++;
129 LOG(2, " count= %2d | channel index= LW_%03d | nu= %.4f cm^-1", n,
130 list_15mu_low[i], l1.nu_lw[list_15mu_low[i]]);
131 }
132 LOG(2, " nu_mean= %.4f cm^-1", nu / n);
133
134 /* 15 micron high channels... */
135 LOG(2, "15 micron high channels:");
136 n = 0;
137 nu = 0;
138 for (int i = 0; i < N15_HIGH; i++) {
139 nu += l1.nu_lw[list_15mu_high[i]];
140 n++;
141 LOG(2, " count= %2d | channel index= LW_%03d | nu= %.4f cm^-1", n,
142 list_15mu_high[i], l1.nu_lw[list_15mu_high[i]]);
143 }
144 LOG(2, " nu_mean= %.4f cm^-1", nu / n);
145
146 /* 8.1 micron channels... */
147 LOG(2, "8.1 micron channels:");
148 n = 0;
149 nu = 0;
150 for (int i = 0; i < N8; i++) {
151 nu += l1.nu_mw[list_8mu[i]];
152 n++;
153 LOG(2, " count= %2d | channel index= MW_%03d | nu= %.4f cm^-1", n,
154 list_8mu[i], l1.nu_mw[list_8mu[i]]);
155 }
156 LOG(2, " nu_mean= %.4f cm^-1", nu / n);
157
158 /* 10.4 micron channels... */
159 LOG(2, "10.4 micron channels:");
160 n = 0;
161 nu = 0;
162 for (int i = 0; i < N10; i++) {
163 nu += l1.nu_lw[list_10mu[i]];
164 n++;
165 LOG(2, " count= %2d | channel index= LW_%03d | nu= %.4f cm^-1", n,
166 list_10mu[i], l1.nu_lw[list_10mu[i]]);
167 }
168 LOG(2, " nu_mean= %.4f cm^-1", nu / n);
169 }
170
171 /* Save geolocation... */
172 pert_4mu->ntrack += L1_NTRACK;
173 if (pert_4mu->ntrack > PERT_NTRACK)
174 ERRMSG("Too many granules!");
175 pert_4mu->nxtrack = L1_NXTRACK;
176 if (pert_4mu->nxtrack > PERT_NXTRACK)
177 ERRMSG("Too many tracks!");
178 pert_4mu->nfov = L1_NFOV;
179 if (pert_4mu->nfov > PERT_NFOV)
180 ERRMSG("Too many field of views!");
181 LOOP_ALL(L1_NTRACK) {
182 pert_4mu->time[track0 + track][xtrack][ifov]
183 = l1.time[track][xtrack] - 220838400.;
184 pert_4mu->lon[track0 + track][xtrack][ifov]
185 = l1.lon[track][xtrack][ifov];
186 pert_4mu->lat[track0 + track][xtrack][ifov]
187 = l1.lat[track][xtrack][ifov];
188 }
189
190 /* Get 8.1 micron brightness temperature... */
191 LOOP_ALL(L1_NTRACK)
192 pert_4mu->dc[track0 + track][xtrack][ifov]
193 = BRIGHT(l1.rad_mw[track][xtrack][ifov][list_8mu[0]] * 1e-3,
194 l1.nu_mw[list_8mu[0]]);
195
196 /* Get 10.4 micron brightness temperature... */
197 LOOP_ALL(L1_NTRACK)
198 pert_15mu_high->dc[track0 + track][xtrack][ifov]
199 = BRIGHT(l1.rad_lw[track][xtrack][ifov][list_10mu[0]] * 1e-3,
200 l1.nu_lw[list_10mu[0]]);
201
202 /* Get 4.3 micron brightness temperature... */
203 LOOP_ALL(L1_NTRACK) {
204 n = 0;
205 nu = rad = 0;
206 for (int i = 0; i < N4; i++)
207 if (gsl_finite(l1.rad_sw[track][xtrack][ifov][list_4mu[i]])) {
208 rad += l1.rad_sw[track][xtrack][ifov][list_4mu[i]] * 1e-3;
209 nu += l1.nu_sw[list_4mu[i]];
210 n++;
211 }
212 if (n > 0.9 * N4)
213 pert_4mu->bt[track0 + track][xtrack][ifov] = BRIGHT(rad / n, nu / n);
214 else
215 pert_4mu->bt[track0 + track][xtrack][ifov] = GSL_NAN;
216 }
217
218 /* Get 15 micron low brightness temperature... */
219 LOOP_ALL(L1_NTRACK) {
220 n = 0;
221 nu = rad = 0;
222 for (int i = 0; i < N15_LOW; i++)
223 if (gsl_finite(l1.rad_lw[track][xtrack][ifov][list_15mu_low[i]])) {
224 rad += l1.rad_lw[track][xtrack][ifov][list_15mu_low[i]] * 1e-3;
225 nu += l1.nu_lw[list_15mu_low[i]];
226 n++;
227 }
228 if (n > 0.9 * N15_LOW)
229 pert_15mu_low->bt[track0 + track][xtrack][ifov] =
230 BRIGHT(rad / n, nu / n);
231 else
232 pert_15mu_low->bt[track0 + track][xtrack][ifov] = GSL_NAN;
233 }
234
235 /* Get 15 micron high brightness temperature... */
236 LOOP_ALL(L1_NTRACK) {
237 n = 0;
238 nu = rad = 0;
239 for (int i = 0; i < N15_HIGH; i++)
240 if (gsl_finite(l1.rad_lw[track][xtrack][ifov][list_15mu_high[i]])) {
241 rad += l1.rad_lw[track][xtrack][ifov][list_15mu_high[i]] * 1e-3;
242 nu += l1.nu_lw[list_15mu_high[i]];
243 n++;
244 }
245 if (n > 0.9 * N15_HIGH)
246 pert_15mu_high->bt[track0 + track][xtrack][ifov] =
247 BRIGHT(rad / n, nu / n);
248 else
249 pert_15mu_high->bt[track0 + track][xtrack][ifov] = GSL_NAN;
250 }
251
252 /* Increment track counter... */
253 track0 += L1_NTRACK;
254 }
255
256 /* ------------------------------------------------------------
257 Outlier filter...
258 ------------------------------------------------------------ */
259
260 /* Write info... */
261 LOG(1, "Outlier filter...");
262
263 /* Define outlier filter... */
264#define OUTLIER_FILTER(BT, DTMED) \
265 do { \
266 for (int track = 0; track < (BT)->ntrack; track++) \
267 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++) { \
268 \
269 size_t nf = 0, outlier_count = 0; \
270 double xf[L1_NFOV]; \
271 \
272 /* Calculate median... */ \
273 for (int ifov = 0; ifov < L1_NFOV; ifov++) { \
274 double v = (BT)->bt[track][xtrack][ifov]; \
275 if (isfinite(v)) xf[nf++] = v; \
276 } \
277 if (nf < 2) continue; \
278 double median = gsl_stats_median(xf, 1, nf); \
279 \
280 /* First pass: count potential outliers... */ \
281 for (int ifov = 0; ifov < L1_NFOV; ifov++) { \
282 double v = (BT)->bt[track][xtrack][ifov]; \
283 if (!isfinite(v)) continue; \
284 if (fabs(v - median) >= (DTMED)) \
285 outlier_count++; \
286 } \
287 \
288 /* Only apply filtering if 1 or 2 outliers exist... */ \
289 if (outlier_count > 0 && outlier_count <= 2) { \
290 for (int ifov = 0; ifov < L1_NFOV; ifov++) { \
291 double v = (BT)->bt[track][xtrack][ifov]; \
292 if (!isfinite(v)) continue; \
293 if (fabs(v - median) >= (DTMED)) { \
294 LOG(2, "outlier: track=%d | xtrack=%d | ifov=%d | " \
295 "bt_" #BT "=%g K | median=%g K", \
296 track, xtrack, ifov, v, median); \
297 (BT)->bt[track][xtrack][ifov] = NAN; \
298 } \
299 } \
300 } \
301 } \
302 } while (0)
303
304 /* Apply filter to each dataset... */
305 OUTLIER_FILTER(pert_4mu, dtmed);
306 OUTLIER_FILTER(pert_15mu_low, dtmed);
307 OUTLIER_FILTER(pert_15mu_high, dtmed);
308
309 /* ------------------------------------------------------------
310 Calculate perturbations...
311 ------------------------------------------------------------ */
312
313 /* Write info... */
314 LOG(1, "Calculate perturbations...");
315
316 /* Loop over scans and field of views... */
317 for (int track = 0; track < pert_4mu->ntrack; track++)
318 for (int ifov = 0; ifov < L1_NFOV; ifov++) {
319
320 /* Skip scan edges... */
321 pert_4mu->dc[track][0][ifov] = GSL_NAN;
322 pert_4mu->dc[track][L1_NXTRACK - 1][ifov] = GSL_NAN;
323 pert_15mu_high->dc[track][0][ifov] = GSL_NAN;
324 pert_15mu_high->dc[track][L1_NXTRACK - 1][ifov] = GSL_NAN;
325 pert_4mu->bt[track][0][ifov] = GSL_NAN;
326 pert_4mu->bt[track][L1_NXTRACK - 1][ifov] = GSL_NAN;
327 pert_15mu_low->bt[track][0][ifov] = GSL_NAN;
328 pert_15mu_low->bt[track][L1_NXTRACK - 1][ifov] = GSL_NAN;
329 pert_15mu_high->bt[track][0][ifov] = GSL_NAN;
330 pert_15mu_high->bt[track][L1_NXTRACK - 1][ifov] = GSL_NAN;
331
332 /* Get 4 micron perturbations... */
333 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
334 x[xtrack] = (double) xtrack;
335 y[xtrack] = pert_4mu->bt[track][xtrack][ifov];
336 }
337 background_poly_help(x, y, L1_NXTRACK, 5);
338 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++)
339 pert_4mu->pt[track][xtrack][ifov] =
340 pert_4mu->bt[track][xtrack][ifov] - y[xtrack];
341
342 /* Get 15 micron low perturbations... */
343 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
344 x[xtrack] = (double) xtrack;
345 y[xtrack] = pert_15mu_low->bt[track][xtrack][ifov];
346 }
347 background_poly_help(x, y, L1_NXTRACK, 7);
348 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++)
349 pert_15mu_low->pt[track][xtrack][ifov] =
350 pert_15mu_low->bt[track][xtrack][ifov] - y[xtrack];
351
352 /* Get 15 micron high perturbations... */
353 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
354 x[xtrack] = (double) xtrack;
355 y[xtrack] = pert_15mu_high->bt[track][xtrack][ifov];
356 }
357 background_poly_help(x, y, L1_NXTRACK, 5);
358 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++)
359 pert_15mu_high->pt[track][xtrack][ifov] =
360 pert_15mu_high->bt[track][xtrack][ifov] - y[xtrack];
361 }
362
363 /* ------------------------------------------------------------
364 Bias correction...
365 ------------------------------------------------------------ */
366
367 /* Check flag... */
368 if (bias) {
369
370 /* Write info... */
371 LOG(1, "Calculate bias correction...");
372
373 /* Radius of Gaussian filter... */
374 const int radius = 25;
375
376 /* Gaussian sigma chosen to yield 50% sensitivity cutoff over ~1000 km... */
377 const double sigma = 11.0;
378 const double inv2sigma2 = 1.0 / (2.0 * sigma * sigma);
379
380 /* Macro for bias correction using Gaussian filter */
381#define APPLY_BIAS_CORR(pert) \
382 do { \
383 LOOP_ALL((pert)->ntrack) { \
384 const int idx = (track * L1_NXTRACK + xtrack) * L1_NFOV + ifov; \
385 help[idx] = 0.0; \
386 double wsum = 0.0; \
387 for (int dtrack2 = -radius; dtrack2 <= radius; ++dtrack2) { \
388 int t2 = track + dtrack2; \
389 if (t2 < 0 || t2 >= (pert)->ntrack) \
390 continue; \
391 if (!isfinite((pert)->pt[t2][xtrack][ifov])) \
392 continue; \
393 const double w = exp(-(double)(dtrack2 * dtrack2) * inv2sigma2); \
394 help[idx] += w * (pert)->pt[t2][xtrack][ifov]; \
395 wsum += w; \
396 } \
397 if (wsum > 0.0) help[idx] /= wsum; \
398 } \
399 LOOP_ALL((pert)->ntrack) { \
400 const int idx = (track * L1_NXTRACK + xtrack) * L1_NFOV + ifov; \
401 (pert)->pt[track][xtrack][ifov] -= help[idx]; \
402 } \
403 } while (0)
404
405 /* Apply bias correction... */
406 APPLY_BIAS_CORR(pert_4mu);
407 APPLY_BIAS_CORR(pert_15mu_low);
408 APPLY_BIAS_CORR(pert_15mu_high);
409 }
410
411 /* ------------------------------------------------------------
412 Calculate variances...
413 ------------------------------------------------------------ */
414
415 /* Write info... */
416 LOG(1, "Calculate variances...");
417
418 /* Loop over footprints... */
419 LOOP_ALL(pert_4mu->ntrack) {
420
421 /* Get geolocation... */
422 geo2cart(0.0, pert_4mu->lon[track][xtrack][ifov],
423 pert_4mu->lat[track][xtrack][ifov], x0);
424
425 /* Init... */
426 int n_4mu = 0, n_15mu_low = 0, n_15mu_high = 0;
427 double mean_4mu = 0, mean_15mu_low = 0, mean_15mu_high = 0;
428 double var_4mu = 0, var_15mu_low = 0, var_15mu_high = 0;
429
430 /* Loop over neighbouring footprints... */
431 for (int track2 = track - dtrack; track2 <= track + dtrack; track2++)
432 for (int xtrack2 = xtrack - dxtrack; xtrack2 <= xtrack + dxtrack;
433 xtrack2++)
434 for (int ifov2 = 0; ifov2 < L1_NFOV; ifov2++)
435 if (track2 >= 0 && track2 < pert_4mu->ntrack
436 && xtrack2 >= 0 && xtrack2 < L1_NXTRACK) {
437
438 /* Check horizontal distance... */
439 geo2cart(0.0, pert_4mu->lon[track2][xtrack2][ifov2],
440 pert_4mu->lat[track2][xtrack2][ifov2], x1);
441 if (DIST2(x0, x1) > var_dh2)
442 continue;
443
444 /* Calculate variance... */
445 if (gsl_finite(pert_4mu->pt[track2][xtrack2][ifov2])) {
446 mean_4mu += pert_4mu->pt[track2][xtrack2][ifov2];
447 var_4mu += gsl_pow_2(pert_4mu->pt[track2][xtrack2][ifov2]);
448 n_4mu++;
449 }
450
451 if (gsl_finite(pert_15mu_low->pt[track2][xtrack2][ifov2])) {
452 mean_15mu_low += pert_15mu_low->pt[track2][xtrack2][ifov2];
453 var_15mu_low +=
454 gsl_pow_2(pert_15mu_low->pt[track2][xtrack2][ifov2]);
455 n_15mu_low++;
456 }
457
458 if (gsl_finite(pert_15mu_high->pt[track2][xtrack2][ifov2])) {
459 mean_15mu_high += pert_15mu_high->pt[track2][xtrack2][ifov2];
460 var_15mu_high +=
461 gsl_pow_2(pert_15mu_high->pt[track2][xtrack2][ifov2]);
462 n_15mu_high++;
463 }
464 }
465
466 /* Save variance data... */
467 if (n_4mu > 0 && xtrack > 0 && xtrack < L1_NXTRACK - 1)
468 pert_4mu->var[track][xtrack][ifov] =
469 var_4mu / n_4mu - gsl_pow_2(mean_4mu / n_4mu);
470 else
471 pert_4mu->var[track][xtrack][ifov] = GSL_NAN;
472
473 if (n_15mu_low > 0 && xtrack > 0 && xtrack < L1_NXTRACK - 1)
474 pert_15mu_low->var[track][xtrack][ifov] =
475 var_15mu_low / n_15mu_low - gsl_pow_2(mean_15mu_low / n_15mu_low);
476 else
477 pert_15mu_low->var[track][xtrack][ifov] = GSL_NAN;
478
479 if (n_15mu_high > 0 && xtrack > 0 && xtrack < L1_NXTRACK - 1)
480 pert_15mu_high->var[track][xtrack][ifov] =
481 var_15mu_high / n_15mu_high - gsl_pow_2(mean_15mu_high / n_15mu_high);
482 else
483 pert_15mu_high->var[track][xtrack][ifov] = GSL_NAN;
484 }
485
486 /* ------------------------------------------------------------
487 Write netCDF file...
488 ------------------------------------------------------------ */
489
490 /* Write info... */
491 LOG(1, "Write perturbation data file: %s", argv[2]);
492
493 /* Create netCDF file... */
494 NC(nc_create(argv[2], NC_CLOBBER, &ncid));
495
496 /* Set dimensions... */
497 NC(nc_def_dim(ncid, "NTRACK", (size_t) pert_4mu->ntrack, &dimid[0]));
498 NC(nc_def_dim(ncid, "NXTRACK", L1_NXTRACK, &dimid[1]));
499 NC(nc_def_dim(ncid, "NFOV", L1_NFOV, &dimid[2]));
500
501 /* Add variables... */
502 NC(nc_def_var(ncid, "time", NC_DOUBLE, 3, dimid, &time_varid));
503 add_att(ncid, time_varid, "s", "time (seconds since 2000-01-01T00:00Z)");
504 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 3, dimid, &lon_varid));
505 add_att(ncid, lon_varid, "deg", "footprint longitude");
506 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 3, dimid, &lat_varid));
507 add_att(ncid, lat_varid, "deg", "footprint latitude");
508
509 NC(nc_def_var(ncid, "bt_8mu", NC_FLOAT, 3, dimid, &bt_8mu_varid));
510 add_att(ncid, bt_8mu_varid, "K", "brightness temperature at 8.1 micron");
511
512 NC(nc_def_var(ncid, "bt_10mu", NC_FLOAT, 3, dimid, &bt_10mu_varid));
513 add_att(ncid, bt_10mu_varid, "K", "brightness temperature at 10.4 micron");
514
515 NC(nc_def_var(ncid, "bt_4mu", NC_FLOAT, 3, dimid, &bt_4mu_varid));
516 add_att(ncid, bt_4mu_varid, "K", "brightness temperature" " at 4.3 micron");
517 NC(nc_def_var(ncid, "bt_4mu_pt", NC_FLOAT, 3, dimid, &bt_4mu_pt_varid));
518 add_att(ncid, bt_4mu_pt_varid, "K", "brightness temperature perturbation"
519 " at 4.3 micron");
520 NC(nc_def_var(ncid, "bt_4mu_var", NC_FLOAT, 3, dimid, &bt_4mu_var_varid));
521 add_att(ncid, bt_4mu_var_varid, "K^2", "brightness temperature variance"
522 " at 4.3 micron");
523
524 NC(nc_def_var(ncid, "bt_15mu_low", NC_FLOAT, 3, dimid, &bt_15mu_low_varid));
525 add_att(ncid, bt_15mu_low_varid, "K", "brightness temperature"
526 " at 15 micron (low altitudes)");
527 NC(nc_def_var(ncid, "bt_15mu_low_pt", NC_FLOAT, 3, dimid,
528 &bt_15mu_low_pt_varid));
529 add_att(ncid, bt_15mu_low_pt_varid, "K",
530 "brightness temperature perturbation"
531 " at 15 micron (low altitudes)");
532 NC(nc_def_var
533 (ncid, "bt_15mu_low_var", NC_FLOAT, 3, dimid, &bt_15mu_low_var_varid));
534 add_att(ncid, bt_15mu_low_var_varid, "K^2",
535 "brightness temperature variance" " at 15 micron (low altitudes)");
536
537 NC(nc_def_var(ncid, "bt_15mu_high", NC_FLOAT, 3, dimid,
538 &bt_15mu_high_varid));
539 add_att(ncid, bt_15mu_high_varid, "K", "brightness temperature"
540 " at 15 micron (high altitudes)");
541 NC(nc_def_var(ncid, "bt_15mu_high_pt", NC_FLOAT, 3, dimid,
542 &bt_15mu_high_pt_varid));
543 add_att(ncid, bt_15mu_high_pt_varid, "K",
544 "brightness temperature perturbation"
545 " at 15 micron (high altitudes)");
546 NC(nc_def_var
547 (ncid, "bt_15mu_high_var", NC_FLOAT, 3, dimid, &bt_15mu_high_var_varid));
548 add_att(ncid, bt_15mu_high_var_varid, "K^2",
549 "brightness temperature variance" " at 15 micron (high altitudes)");
550
551 /* Leave define mode... */
552 NC(nc_enddef(ncid));
553
554 /* Write data... */
555 WRITE_FIELD(pert_4mu->time[track][xtrack][ifov], time_varid);
556 WRITE_FIELD(pert_4mu->lon[track][xtrack][ifov], lon_varid);
557 WRITE_FIELD(pert_4mu->lat[track][xtrack][ifov], lat_varid);
558 WRITE_FIELD(pert_4mu->dc[track][xtrack][ifov], bt_8mu_varid);
559 WRITE_FIELD(pert_15mu_high->dc[track][xtrack][ifov], bt_10mu_varid);
560 WRITE_FIELD(pert_4mu->bt[track][xtrack][ifov], bt_4mu_varid);
561 WRITE_FIELD(pert_4mu->pt[track][xtrack][ifov], bt_4mu_pt_varid);
562 WRITE_FIELD(pert_4mu->var[track][xtrack][ifov], bt_4mu_var_varid);
563 WRITE_FIELD(pert_15mu_low->bt[track][xtrack][ifov], bt_15mu_low_varid);
564 WRITE_FIELD(pert_15mu_low->pt[track][xtrack][ifov], bt_15mu_low_pt_varid);
565 WRITE_FIELD(pert_15mu_low->var[track][xtrack][ifov], bt_15mu_low_var_varid);
566 WRITE_FIELD(pert_15mu_high->bt[track][xtrack][ifov], bt_15mu_high_varid);
567 WRITE_FIELD(pert_15mu_high->pt[track][xtrack][ifov], bt_15mu_high_pt_varid);
568 WRITE_FIELD(pert_15mu_high->var[track][xtrack][ifov],
569 bt_15mu_high_var_varid);
570
571 /* Close file... */
572 NC(nc_close(ncid));
573
574 /* Free... */
575 free(pert_4mu);
576 free(pert_15mu_low);
577 free(pert_15mu_high);
578
579 return EXIT_SUCCESS;
580}
scan_ctl
double scan_ctl(int argc, char *argv[], const char *varname, const int arridx, const char *defvalue, char *value)
Scan control file or command-line arguments for a configuration variable.
Definition: jurassic.c:6267
geo2cart
void geo2cart(const double z, const double lon, const double lat, double *x)
Converts geographic coordinates (longitude, latitude, altitude) to Cartesian coordinates.
Definition: jurassic.c:3853
BRIGHT
#define BRIGHT(rad, nu)
Compute brightness temperature from radiance.
Definition: jurassic.h:412
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: jurassic.h:948
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for an array.
Definition: jurassic.h:382
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: jurassic.h:878
DIST2
#define DIST2(a, b)
Compute squared distance between two 3D vectors.
Definition: jurassic.h:464
background_poly_help
void background_poly_help(const double *xx, double *yy, const int n, const int dim)
Get background based on polynomial fits.
Definition: libcris.c:47
read_cris_l1
int read_cris_l1(char *filename, cris_l1_t *l1, int apo)
Read CrIS Level-1 data.
Definition: libcris.c:1051
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
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libcris.h:61
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libcris.h:42
L1_NTRACK
#define L1_NTRACK
Along-track size of CrIS radiance granule.
Definition: libcris.h:21
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libcris.h:39
PERT_NFOV
#define PERT_NFOV
Number of field of views of perturbation data.
Definition: libcris.h:45
OUTLIER_FILTER
#define OUTLIER_FILTER(BT, DTMED)
N15_LOW
#define N15_LOW
Number of 15 micron low channels:
Definition: perturbation.c:11
N10
#define N10
Number of 10.4 micron channels:
Definition: perturbation.c:20
WRITE_FIELD
#define WRITE_FIELD(expr, varid)
Definition: perturbation.c:33
APPLY_BIAS_CORR
#define APPLY_BIAS_CORR(pert)
N4
#define N4
Number of 4.3 micron channels:
Definition: perturbation.c:8
N8
#define N8
Number of 8.1 micron channels:
Definition: perturbation.c:17
N15_HIGH
#define N15_HIGH
Number of 15 micron high channels:
Definition: perturbation.c:14
cris_l1_t
CrIS Level-1 data.
Definition: libcris.h:72
cris_l1_t::nu_sw
double nu_sw[L1_NCHAN_SW]
Shortwave channel frequencies [cm^-1].
Definition: libcris.h:117
cris_l1_t::rad_mw
float rad_mw[L1_NTRACK][L1_NXTRACK][L1_NFOV][L1_NCHAN_MW]
Midwave radiance [W/(m^2 sr cm^-1)].
Definition: libcris.h:108
cris_l1_t::nu_mw
double nu_mw[L1_NCHAN_MW]
Midwave channel frequencies [cm^-1].
Definition: libcris.h:105
cris_l1_t::rad_sw
float rad_sw[L1_NTRACK][L1_NXTRACK][L1_NFOV][L1_NCHAN_SW]
Shortwave radiance [W/(m^2 sr cm^-1)].
Definition: libcris.h:120
cris_l1_t::lon
double lon[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Footprint longitude [deg].
Definition: libcris.h:78
cris_l1_t::rad_lw
float rad_lw[L1_NTRACK][L1_NXTRACK][L1_NFOV][L1_NCHAN_LW]
Longwave radiance [W/(m^2 sr cm^-1)].
Definition: libcris.h:96
cris_l1_t::time
double time[L1_NTRACK][L1_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libcris.h:75
cris_l1_t::nu_lw
double nu_lw[L1_NCHAN_LW]
Longwave channel frequencies [cm^-1].
Definition: libcris.h:93
cris_l1_t::lat
double lat[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Footprint latitude [deg].
Definition: libcris.h:81
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
Here is the call graph for this function:
Generated by doxygen 1.9.4