IASI Code Collection
Macros | Functions
perturbation.c File Reference

Extract perturbation data. More...

#include "libiasi.h"

Go to the source code of this file.

Macros

#define N4   152
 
#define N15_LOW   21
 
#define N15_HIGH   2
 

Functions

void addatt (int ncid, int varid, const char *unit, const char *long_name)
 
int main (int argc, char *argv[])
 

Detailed Description

Extract perturbation data.

Definition in file perturbation.c.

Macro Definition Documentation

◆ N4

#define N4   152

Definition at line 33 of file perturbation.c.

◆ N15_LOW

#define N15_LOW   21

Definition at line 36 of file perturbation.c.

◆ N15_HIGH

#define N15_HIGH   2

Definition at line 39 of file perturbation.c.

Function Documentation

◆ addatt()

void addatt ( int  ncid,
int  varid,
const char *  unit,
const char *  long_name 
)

Definition at line 434 of file perturbation.c.

438 {
439
440 /* Set long name... */
441 NC(nc_put_att_text(ncid, varid, "long_name", strlen(long_name), long_name));
442
443 /* Set units... */
444 NC(nc_put_att_text(ncid, varid, "units", strlen(unit), unit));
445}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libiasi.h:155

◆ main()

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

Definition at line 56 of file perturbation.c.

58 {
59
60 static iasi_rad_t *iasi_rad;
61
62 static pert_t *pert_4mu, *pert_15mu_low, *pert_15mu_high;
63
64 static wave_t wave;
65
66 static double numean, radmean, var_dh = 100.;
67
68 const int list_4mu[N4]
69 = { 6711, 6712, 6713, 6714, 6715, 6716, 6717, 6718, 6719, 6720,
70 6721, 6722, 6723, 6724, 6725, 6726, 6727, 6728, 6729, 6730, 6731,
71 6732, 6733, 6734, 6735, 6736, 6737, 6738, 6739, 6740, 6741, 6742,
72 6743, 6744, 6745, 6746, 6747, 6748, 6749, 6750, 6751, 6752, 6753,
73 6754, 6755, 6756, 6757, 6758, 6759, 6760, 6761, 6762, 6763, 6764,
74 6765, 6766, 6767, 6768, 6769, 6770, 6771, 6772, 6773, 6774, 6775,
75 6776, 6777, 6778, 6779, 6780, 6781, 6782, 6783, 6784, 6785, 6786,
76 6787, 6788, 6789, 6790, 6791, 6792, 6793, 6794, 6795, 6796, 6797,
77 6798, 6799, 6800, 6801, 6802, 6803, 6804, 6830, 6831, 6832, 6833,
78 6834, 6835, 6836, 6837, 6838, 6839, 6840, 6841, 6842, 6843, 6844,
79 6845, 6846, 6847, 6848, 6849, 6850, 6851, 6852, 6853, 6854, 6855,
80 6856, 6857, 6858, 6859, 6860, 6861, 6862, 6863, 6864, 6865, 6866,
81 6867, 6868, 6869, 6870, 6871, 6872, 6873, 6874, 6875, 6876, 6877,
82 6878, 6879, 6880, 6881, 6882, 6883, 6884, 6885, 6886, 6887
83 };
84
85 const int list_15mu_low[N15_LOW]
86 = { 22, 28, 34, 40, 46, 52, 58, 72, 100, 105, 112, 118, 119,
87 124, 125, 130, 131, 136, 137, 143, 144
88 };
89
90 const int list_15mu_high[N15_HIGH]
91 = { 91, 92 };
92
93 const int cloud_chan = 2364;
94
95 static int ix, iy, dimid[2], i, n, ncid, track, track0, xtrack,
96 time_varid, lon_varid, lat_varid, bt_4mu_varid, bt_4mu_pt_varid,
97 bt_4mu_var_varid, bt_8mu_varid, bt_15mu_low_varid, bt_15mu_low_pt_varid,
98 bt_15mu_low_var_varid, bt_15mu_high_varid, bt_15mu_high_pt_varid,
99 bt_15mu_high_var_varid, iarg, format, init;
100
101 static size_t start[2], count[2];
102
103 /* Check arguments... */
104 if (argc < 4)
105 ERRMSG("Give parameters: <ctl> <out.nc> <l1b_file1> [<l1b_file2> ...]");
106
107 /* Read control parameters... */
108 format = (int) scan_ctl(argc, argv, "FORMAT", -1, "1", NULL);
109
110 /* Allocate... */
111 ALLOC(iasi_rad, iasi_rad_t, 1);
112 ALLOC(pert_4mu, pert_t, 1);
113 ALLOC(pert_15mu_low, pert_t, 1);
114 ALLOC(pert_15mu_high, pert_t, 1);
115
116 /* ------------------------------------------------------------
117 Read HDF files...
118 ------------------------------------------------------------ */
119
120 /* Loop over IASI files... */
121 for (iarg = 3; iarg < argc; iarg++) {
122
123 /* Read IASI data... */
124 printf("Read IASI Level-1C data file: %s\n", argv[iarg]);
125 iasi_read(format, argv[iarg], iasi_rad);
126
127 /* Write info... */
128 if (!init) {
129 init = 1;
130 LOG(2, "4 micron channels:");
131 for (i = 0; i < N4; i++)
132 LOG(2, " channel[%4d]= %4d | freq[%4d]= %7.2f cm^-1", i, list_4mu[i],
133 i, iasi_rad->freq[list_4mu[i]]);
134 LOG(2, "15 micron low channels:");
135 for (i = 0; i < N15_LOW; i++)
136 LOG(2, " channel[%4d]= %4d | freq[%4d]= %7.2f cm^-1", i,
137 list_15mu_low[i], i, iasi_rad->freq[list_15mu_low[i]]);
138 LOG(2, "15 micron high channels:");
139 for (i = 0; i < N15_HIGH; i++)
140 LOG(2, " channel[%4d]= %4d | freq[%4d]= %7.2f cm^-1", i,
141 list_15mu_high[i], i, iasi_rad->freq[list_15mu_high[i]]);
142 LOG(2, "cloud channel:");
143 LOG(2, " channel[%4d]= %4d | freq[%4d]= %7.2f cm^-1", 0, cloud_chan, 0,
144 iasi_rad->freq[cloud_chan]);
145 }
146
147 /* Save geolocation... */
148 pert_4mu->ntrack += iasi_rad->ntrack;
149 if (pert_4mu->ntrack > PERT_NTRACK)
150 ERRMSG("Too many granules!");
151 pert_4mu->nxtrack = L1_NXTRACK;
152 if (pert_4mu->nxtrack > PERT_NXTRACK)
153 ERRMSG("Too many tracks!");
154 for (track = 0; track < iasi_rad->ntrack; track++)
155 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
156 pert_4mu->time[track0 + track][xtrack]
157 = iasi_rad->Time[track][xtrack];
158 pert_4mu->lon[track0 + track][xtrack]
159 = iasi_rad->Longitude[track][xtrack];
160 pert_4mu->lat[track0 + track][xtrack]
161 = iasi_rad->Latitude[track][xtrack];
162 }
163
164 pert_15mu_low->ntrack += iasi_rad->ntrack;
165 if (pert_15mu_low->ntrack > PERT_NTRACK)
166 ERRMSG("Too many granules!");
167 pert_15mu_low->nxtrack = L1_NXTRACK;
168 if (pert_15mu_low->nxtrack > PERT_NXTRACK)
169 ERRMSG("Too many tracks!");
170 for (track = 0; track < iasi_rad->ntrack; track++)
171 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
172 pert_15mu_low->time[track0 + track][xtrack]
173 = iasi_rad->Time[track][xtrack];
174 pert_15mu_low->lon[track0 + track][xtrack]
175 = iasi_rad->Longitude[track][xtrack];
176 pert_15mu_low->lat[track0 + track][xtrack]
177 = iasi_rad->Latitude[track][xtrack];
178 }
179
180 pert_15mu_high->ntrack += iasi_rad->ntrack;
181 if (pert_15mu_high->ntrack > PERT_NTRACK)
182 ERRMSG("Too many granules!");
183 pert_15mu_high->nxtrack = L1_NXTRACK;
184 if (pert_15mu_high->nxtrack > PERT_NXTRACK)
185 ERRMSG("Too many tracks!");
186 for (track = 0; track < iasi_rad->ntrack; track++)
187 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
188 pert_15mu_high->time[track0 + track][xtrack]
189 = iasi_rad->Time[track][xtrack];
190 pert_15mu_high->lon[track0 + track][xtrack]
191 = iasi_rad->Longitude[track][xtrack];
192 pert_15mu_high->lat[track0 + track][xtrack]
193 = iasi_rad->Latitude[track][xtrack];
194 }
195
196 /* Get 8.1 micron brightness temperature... */
197 for (track = 0; track < iasi_rad->ntrack; track++)
198 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++)
199 pert_4mu->dc[track0 + track][xtrack]
200 = BRIGHT(iasi_rad->Rad[track][xtrack][cloud_chan],
201 iasi_rad->freq[cloud_chan]);
202
203 /* Get 4.3 micron brightness temperature... */
204 for (track = 0; track < iasi_rad->ntrack; track++)
205 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
206 n = 0;
207 numean = radmean = 0;
208 for (i = 0; i < N4; i++)
209 if (gsl_finite(iasi_rad->Rad[track][xtrack][list_4mu[i]])) {
210 radmean += iasi_rad->Rad[track][xtrack][list_4mu[i]];
211 numean += iasi_rad->freq[list_4mu[i]];
212 n++;
213 }
214 if (n > 0.9 * N4)
215 pert_4mu->bt[track0 + track][xtrack]
216 = BRIGHT(radmean / n, numean / n);
217 else
218 pert_4mu->bt[track0 + track][xtrack] = GSL_NAN;
219 }
220
221 /* Get 15 micron brightness temperature (low altitudes)... */
222 for (track = 0; track < iasi_rad->ntrack; track++)
223 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
224 n = 0;
225 numean = radmean = 0;
226 for (i = 0; i < N15_LOW; i++)
227 if (gsl_finite(iasi_rad->Rad[track][xtrack][list_15mu_low[i]])) {
228 radmean += iasi_rad->Rad[track][xtrack][list_15mu_low[i]];
229 numean += iasi_rad->freq[list_15mu_low[i]];
230 n++;
231 }
232 if (n > 0.9 * N15_LOW)
233 pert_15mu_low->bt[track0 + track][xtrack]
234 = BRIGHT(radmean / n, numean / n);
235 else
236 pert_15mu_low->bt[track0 + track][xtrack] = GSL_NAN;
237 }
238
239 /* Get 15 micron brightness temperature (high altitudes)... */
240 for (track = 0; track < iasi_rad->ntrack; track++)
241 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
242 n = 0;
243 numean = radmean = 0;
244 for (i = 0; i < N15_HIGH; i++)
245 if (gsl_finite(iasi_rad->Rad[track][xtrack][list_15mu_high[i]])) {
246 radmean += iasi_rad->Rad[track][xtrack][list_15mu_high[i]];
247 numean += iasi_rad->freq[list_15mu_high[i]];
248 n++;
249 }
250 if (n > 0.9 * N15_HIGH)
251 pert_15mu_high->bt[track0 + track][xtrack]
252 = BRIGHT(radmean / n, numean / n);
253 else
254 pert_15mu_high->bt[track0 + track][xtrack] = GSL_NAN;
255 }
256
257 /* Increment track counter... */
258 track0 += iasi_rad->ntrack;
259 }
260
261 /* ------------------------------------------------------------
262 Calculate perturbations and variances...
263 ------------------------------------------------------------ */
264
265 /* Convert to wave analysis struct... */
266 pert2wave(pert_4mu, &wave,
267 0, pert_4mu->ntrack - 1, 0, pert_4mu->nxtrack - 1);
268
269 /* Estimate background... */
270 background_poly(&wave, 5, 0);
271
272 /* Compute variance... */
273 variance(&wave, var_dh);
274
275 /* Copy data... */
276 for (ix = 0; ix < wave.nx; ix++)
277 for (iy = 0; iy < wave.ny; iy++) {
278 pert_4mu->pt[iy][ix] = wave.pt[ix][iy];
279 pert_4mu->var[iy][ix] = wave.var[ix][iy];
280 }
281
282 /* Convert to wave analysis struct... */
283 pert2wave(pert_15mu_low, &wave,
284 0, pert_15mu_low->ntrack - 1, 0, pert_15mu_low->nxtrack - 1);
285
286 /* Estimate background... */
287 background_poly(&wave, 5, 0);
288
289 /* Compute variance... */
290 variance(&wave, var_dh);
291
292 /* Copy data... */
293 for (ix = 0; ix < wave.nx; ix++)
294 for (iy = 0; iy < wave.ny; iy++) {
295 pert_15mu_low->pt[iy][ix] = wave.pt[ix][iy];
296 pert_15mu_low->var[iy][ix] = wave.var[ix][iy];
297 }
298
299 /* Convert to wave analysis struct... */
300 pert2wave(pert_15mu_high, &wave,
301 0, pert_15mu_high->ntrack - 1, 0, pert_15mu_high->nxtrack - 1);
302
303 /* Estimate background... */
304 background_poly(&wave, 5, 0);
305
306 /* Compute variance... */
307 variance(&wave, var_dh);
308
309 /* Copy data... */
310 for (ix = 0; ix < wave.nx; ix++)
311 for (iy = 0; iy < wave.ny; iy++) {
312 pert_15mu_high->pt[iy][ix] = wave.pt[ix][iy];
313 pert_15mu_high->var[iy][ix] = wave.var[ix][iy];
314 }
315
316 /* ------------------------------------------------------------
317 Write to netCDF file...
318 ------------------------------------------------------------ */
319
320 /* Create netCDF file... */
321 printf("Write perturbation data file: %s\n", argv[2]);
322 NC(nc_create(argv[2], NC_CLOBBER, &ncid));
323
324 /* Set dimensions... */
325 NC(nc_def_dim(ncid, "NTRACK", NC_UNLIMITED, &dimid[0]));
326 NC(nc_def_dim(ncid, "NXTRACK", L1_NXTRACK, &dimid[1]));
327
328 /* Add variables... */
329 NC(nc_def_var(ncid, "time", NC_DOUBLE, 2, dimid, &time_varid));
330 addatt(ncid, time_varid, "s", "time (seconds since 2000-01-01T00:00Z)");
331 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 2, dimid, &lon_varid));
332 addatt(ncid, lon_varid, "deg", "footprint longitude");
333 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 2, dimid, &lat_varid));
334 addatt(ncid, lat_varid, "deg", "footprint latitude");
335
336 NC(nc_def_var(ncid, "bt_8mu", NC_FLOAT, 2, dimid, &bt_8mu_varid));
337 addatt(ncid, bt_8mu_varid, "K", "brightness temperature at 8.1 micron");
338
339 NC(nc_def_var(ncid, "bt_4mu", NC_FLOAT, 2, dimid, &bt_4mu_varid));
340 addatt(ncid, bt_4mu_varid, "K", "brightness temperature" " at 4.3 micron");
341 NC(nc_def_var(ncid, "bt_4mu_pt", NC_FLOAT, 2, dimid, &bt_4mu_pt_varid));
342 addatt(ncid, bt_4mu_pt_varid, "K", "brightness temperature perturbation"
343 " at 4.3 micron");
344 NC(nc_def_var(ncid, "bt_4mu_var", NC_FLOAT, 2, dimid, &bt_4mu_var_varid));
345 addatt(ncid, bt_4mu_var_varid, "K^2", "brightness temperature variance"
346 " at 4.3 micron");
347
348 NC(nc_def_var(ncid, "bt_15mu_low", NC_FLOAT, 2, dimid, &bt_15mu_low_varid));
349 addatt(ncid, bt_15mu_low_varid, "K", "brightness temperature"
350 " at 15 micron (low altitudes)");
351 NC(nc_def_var(ncid, "bt_15mu_low_pt", NC_FLOAT, 2, dimid,
352 &bt_15mu_low_pt_varid));
353 addatt(ncid, bt_15mu_low_pt_varid, "K",
354 "brightness temperature perturbation"
355 " at 15 micron (low altitudes)");
356 NC(nc_def_var
357 (ncid, "bt_15mu_low_var", NC_FLOAT, 2, dimid, &bt_15mu_low_var_varid));
358 addatt(ncid, bt_15mu_low_var_varid, "K^2",
359 "brightness temperature variance" " at 15 micron (low altitudes)");
360
361 NC(nc_def_var(ncid, "bt_15mu_high", NC_FLOAT, 2, dimid,
362 &bt_15mu_high_varid));
363 addatt(ncid, bt_15mu_high_varid, "K", "brightness temperature"
364 " at 15 micron (high altitudes)");
365 NC(nc_def_var(ncid, "bt_15mu_high_pt", NC_FLOAT, 2, dimid,
366 &bt_15mu_high_pt_varid));
367 addatt(ncid, bt_15mu_high_pt_varid, "K",
368 "brightness temperature perturbation"
369 " at 15 micron (high altitudes)");
370 NC(nc_def_var
371 (ncid, "bt_15mu_high_var", NC_FLOAT, 2, dimid, &bt_15mu_high_var_varid));
372 addatt(ncid, bt_15mu_high_var_varid, "K^2",
373 "brightness temperature variance" " at 15 micron (high altitudes)");
374
375 /* Leave define mode... */
376 NC(nc_enddef(ncid));
377
378 /* Loop over tracks... */
379 for (track = 0; track < pert_4mu->ntrack; track++) {
380
381 /* Set array sizes... */
382 start[0] = (size_t) track;
383 start[1] = 0;
384 count[0] = 1;
385 count[1] = (size_t) pert_4mu->nxtrack;
386
387 /* Write data... */
388 NC(nc_put_vara_double(ncid, time_varid, start, count,
389 pert_4mu->time[track]));
390 NC(nc_put_vara_double(ncid, lon_varid, start, count,
391 pert_4mu->lon[track]));
392 NC(nc_put_vara_double(ncid, lat_varid, start, count,
393 pert_4mu->lat[track]));
394
395 NC(nc_put_vara_double(ncid, bt_8mu_varid, start, count,
396 pert_4mu->dc[track]));
397
398 NC(nc_put_vara_double(ncid, bt_4mu_varid, start, count,
399 pert_4mu->bt[track]));
400 NC(nc_put_vara_double(ncid, bt_4mu_pt_varid, start, count,
401 pert_4mu->pt[track]));
402 NC(nc_put_vara_double(ncid, bt_4mu_var_varid, start, count,
403 pert_4mu->var[track]));
404
405 NC(nc_put_vara_double(ncid, bt_15mu_low_varid, start, count,
406 pert_15mu_low->bt[track]));
407 NC(nc_put_vara_double(ncid, bt_15mu_low_pt_varid, start, count,
408 pert_15mu_low->pt[track]));
409 NC(nc_put_vara_double(ncid, bt_15mu_low_var_varid, start, count,
410 pert_15mu_low->var[track]));
411
412 NC(nc_put_vara_double(ncid, bt_15mu_high_varid, start, count,
413 pert_15mu_high->bt[track]));
414 NC(nc_put_vara_double(ncid, bt_15mu_high_pt_varid, start, count,
415 pert_15mu_high->pt[track]));
416 NC(nc_put_vara_double(ncid, bt_15mu_high_var_varid, start, count,
417 pert_15mu_high->var[track]));
418 }
419
420 /* Close file... */
421 NC(nc_close(ncid));
422
423 /* Free... */
424 free(iasi_rad);
425 free(pert_4mu);
426 free(pert_15mu_low);
427 free(pert_15mu_high);
428
429 return EXIT_SUCCESS;
430}
iasi_read
void iasi_read(int format, char *filename, iasi_rad_t *iasi_rad)
Read IASI Level-1 data.
Definition: libiasi.c:297
background_poly
void background_poly(wave_t *wave, int dim_x, int dim_y)
Get background based on polynomial fits.
Definition: libiasi.c:114
variance
void variance(wave_t *wave, double dh)
Compute local variance.
Definition: libiasi.c:1081
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: libiasi.c:885
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libiasi.h:135
L1_NXTRACK
#define L1_NXTRACK
Across-track size of IASI radiance granule (don't change).
Definition: libiasi.h:102
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libiasi.h:132
N15_LOW
#define N15_LOW
Definition: perturbation.c:36
N4
#define N4
Definition: perturbation.c:33
addatt
void addatt(int ncid, int varid, const char *unit, const char *long_name)
Definition: perturbation.c:434
N15_HIGH
#define N15_HIGH
Definition: perturbation.c:39
iasi_rad_t
IASI converted Level-1 radiation data.
Definition: libiasi.h:288
iasi_rad_t::Longitude
double Longitude[L1_NTRACK][L1_NXTRACK]
Longitude of the sounder pixel.
Definition: libiasi.h:300
iasi_rad_t::Latitude
double Latitude[L1_NTRACK][L1_NXTRACK]
Latitude of the sounder pixel.
Definition: libiasi.h:303
iasi_rad_t::Time
double Time[L1_NTRACK][L1_NXTRACK]
Seconds since 2000-01-01 for each sounder pixel.
Definition: libiasi.h:297
iasi_rad_t::freq
double freq[IASI_L1_NCHAN]
channel wavenumber [cm^-1]
Definition: libiasi.h:294
iasi_rad_t::ntrack
int ntrack
Number of along-track samples.
Definition: libiasi.h:291
iasi_rad_t::Rad
float Rad[L1_NTRACK][L1_NXTRACK][IASI_L1_NCHAN]
Radiance [W/(m^2 sr cm^-1)].
Definition: libiasi.h:306
pert_t
Perturbation data.
Definition: libiasi.h:224
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libiasi.h:251
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libiasi.h:248
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:233
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libiasi.h:242
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libiasi.h:227
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libiasi.h:230
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK]
Longitude [deg].
Definition: libiasi.h:236
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (4 or 15 micron) [K].
Definition: libiasi.h:245
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK]
Latitude [deg].
Definition: libiasi.h:239
wave_t
Wave analysis data.
Definition: libiasi.h:320
wave_t::nx
int nx
Number of across-track values.
Definition: libiasi.h:323
wave_t::ny
int ny
Number of along-track values.
Definition: libiasi.h:326
wave_t::var
double var[WX][WY]
Variance [K].
Definition: libiasi.h:356
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libiasi.h:353
Here is the call graph for this function:
Generated by doxygen 1.9.4