AIRS Code Collection
perturbation.c
Go to the documentation of this file.
1/*
2 This file is part of the AIRS Code Collection.
3
4 the AIRS Code Collections is free software: you can redistribute it
5 and/or modify it under the terms of the GNU General Public License
6 as published by the Free Software Foundation, either version 3 of
7 the License, or (at your option) any later version.
8
9 The AIRS Code Collection is distributed in the hope that it will be
10 useful, but WITHOUT ANY WARRANTY; without even the implied warranty
11 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with the AIRS Code Collection. If not, see
16 <http://www.gnu.org/licenses/>.
17
18 Copyright (C) 2019-2025 Forschungszentrum Juelich GmbH
19*/
20
26#include "libairs.h"
27
28/* ------------------------------------------------------------
29 Constants...
30 ------------------------------------------------------------ */
31
32/* Number of 4 micron channels: */
33#define N4 42
34
35/* Number of 15 micron channels (low altitudes): */
36#define N15_LOW 21
37
38/* Number of 15 micron channels (high altitudes): */
39#define N15_HIGH 2
40
41/* ------------------------------------------------------------
42 Main...
43 ------------------------------------------------------------ */
44
45int main(
46 int argc,
47 char *argv[]) {
48
49 static airs_rad_gran_t airs_rad_gran;
50
51 static pert_t *pert_4mu, *pert_15mu_low, *pert_15mu_high;
52
53 static wave_t wave;
54
55 static double var_dh = 100.;
56
57 static int list_4mu[N4]
58 = { 2039, 2040, 2041, 2042, 2043, 2044, 2045, 2046, 2047, 2048,
59 2049, 2050, 2051, 2052, 2053, 2054, 2055, 2056, 2057, 2058,
60 2059, 2060, 2061, 2062, 2063, 2064, 2071, 2072, 2073, 2074,
61 2075, 2076, 2077, 2078, 2079, 2080, 2081, 2082, 2083, 2084,
62 2085, 2086
63 };
64
65 static int list_15mu_low[N15_LOW]
66 = { 4, 10, 16, 22, 29, 35, 41, 55, 83, 88, 94,
67 100, 101, 106, 107, 112, 113, 118, 119, 124, 125
68 };
69
70 static int list_15mu_high[N15_HIGH]
71 = { 74, 75 };
72
73 static int dimid[2], n, ncid, track0, time_varid, lon_varid, lat_varid,
74 bt_4mu_varid, bt_4mu_pt_varid, bt_4mu_var_varid, bt_8mu_varid,
75 bt_15mu_low_varid, bt_15mu_low_pt_varid, bt_15mu_low_var_varid,
76 bt_15mu_high_varid, bt_15mu_high_pt_varid, bt_15mu_high_var_varid;
77
78 static size_t start[2], count[2];
79
80 /* Check arguments... */
81 if (argc < 3)
82 ERRMSG("Give parameters: <out.nc> <l1b_file1> [<l1b_file2> ...]");
83
84 /* Allocate... */
85 ALLOC(pert_4mu, pert_t, 1);
86 ALLOC(pert_15mu_low, pert_t, 1);
87 ALLOC(pert_15mu_high, pert_t, 1);
88
89 /* ------------------------------------------------------------
90 Read HDF files...
91 ------------------------------------------------------------ */
92
93 /* Loop over HDF files... */
94 for (int iarg = 2; iarg < argc; iarg++) {
95
96 /* Read AIRS data... */
97 printf("Read AIRS Level-1B data file: %s\n", argv[iarg]);
98 airs_rad_rdr(argv[iarg], &airs_rad_gran);
99
100 /* Flag bad observations... */
101 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
102 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
103 for (int i = 0; i < AIRS_RAD_CHANNEL; i++)
104 if ((airs_rad_gran.state[track][xtrack] != 0)
105 || (airs_rad_gran.ExcludedChans[i] > 2)
106 || (airs_rad_gran.CalChanSummary[i] & 8)
107 || (airs_rad_gran.CalChanSummary[i] & (32 + 64))
108 || (airs_rad_gran.CalFlag[track][i] & 16)
109 || (airs_rad_gran.Longitude[track][xtrack] < -180)
110 || (airs_rad_gran.Longitude[track][xtrack] > 180)
111 || (airs_rad_gran.Latitude[track][xtrack] < -90)
112 || (airs_rad_gran.Latitude[track][xtrack] > 90))
113 airs_rad_gran.radiances[track][xtrack][i] = GSL_NAN;
114 else
115 airs_rad_gran.radiances[track][xtrack][i] *= 0.001f;
116
117 /* Save geolocation... */
118 pert_4mu->ntrack += AIRS_RAD_GEOTRACK;
119 if (pert_4mu->ntrack > PERT_NTRACK)
120 ERRMSG("Too many granules!");
121 pert_4mu->nxtrack = AIRS_RAD_GEOXTRACK;
122 if (pert_4mu->nxtrack > PERT_NXTRACK)
123 ERRMSG("Too many tracks!");
124 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
125 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
126 pert_4mu->time[track0 + track][xtrack]
127 = airs_rad_gran.Time[track][xtrack] - 220838400.;
128 pert_4mu->lon[track0 + track][xtrack]
129 = airs_rad_gran.Longitude[track][xtrack];
130 pert_4mu->lat[track0 + track][xtrack]
131 = airs_rad_gran.Latitude[track][xtrack];
132 }
133
134 pert_15mu_low->ntrack += AIRS_RAD_GEOTRACK;
135 if (pert_15mu_low->ntrack > PERT_NTRACK)
136 ERRMSG("Too many granules!");
137 pert_15mu_low->nxtrack = AIRS_RAD_GEOXTRACK;
138 if (pert_15mu_low->nxtrack > PERT_NXTRACK)
139 ERRMSG("Too many tracks!");
140 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
141 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
142 pert_15mu_low->time[track0 + track][xtrack]
143 = airs_rad_gran.Time[track][xtrack] - 220838400.;
144 pert_15mu_low->lon[track0 + track][xtrack]
145 = airs_rad_gran.Longitude[track][xtrack];
146 pert_15mu_low->lat[track0 + track][xtrack]
147 = airs_rad_gran.Latitude[track][xtrack];
148 }
149
150 pert_15mu_high->ntrack += AIRS_RAD_GEOTRACK;
151 if (pert_15mu_high->ntrack > PERT_NTRACK)
152 ERRMSG("Too many granules!");
153 pert_15mu_high->nxtrack = AIRS_RAD_GEOXTRACK;
154 if (pert_15mu_high->nxtrack > PERT_NXTRACK)
155 ERRMSG("Too many tracks!");
156 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
157 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
158 pert_15mu_high->time[track0 + track][xtrack]
159 = airs_rad_gran.Time[track][xtrack] - 220838400.;
160 pert_15mu_high->lon[track0 + track][xtrack]
161 = airs_rad_gran.Longitude[track][xtrack];
162 pert_15mu_high->lat[track0 + track][xtrack]
163 = airs_rad_gran.Latitude[track][xtrack];
164 }
165
166 /* Get 8.1 micron brightness temperature... */
167 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
168 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
169 pert_4mu->dc[track0 + track][xtrack]
170 = BRIGHT(airs_rad_gran.radiances[track][xtrack][1290],
171 airs_rad_gran.nominal_freq[1290]);
172
173 /* Get 4.3 micron brightness temperature... */
174 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
175 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
176 n = 0;
177 for (int i = 0; i < N4; i++)
178 if (gsl_finite(airs_rad_gran.radiances[track][xtrack][list_4mu[i]])) {
179 pert_4mu->bt[track0 + track][xtrack]
180 += BRIGHT(airs_rad_gran.radiances[track][xtrack][list_4mu[i]],
181 airs_rad_gran.nominal_freq[list_4mu[i]]);
182 n++;
183 }
184 if (n > 0.9 * N4)
185 pert_4mu->bt[track0 + track][xtrack] /= n;
186 else
187 pert_4mu->bt[track0 + track][xtrack] = GSL_NAN;
188 }
189
190 /* Get 15 micron brightness temperature (low altitudes)... */
191 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
192 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
193 n = 0;
194 for (int i = 0; i < N15_LOW; i++)
195 if (gsl_finite(airs_rad_gran.radiances
196 [track][xtrack][list_15mu_low[i]])) {
197 pert_15mu_low->bt[track0 + track][xtrack]
198 +=
199 BRIGHT(airs_rad_gran.radiances[track][xtrack][list_15mu_low[i]],
200 airs_rad_gran.nominal_freq[list_15mu_low[i]]);
201 n++;
202 }
203 if (n > 0.9 * N15_LOW)
204 pert_15mu_low->bt[track0 + track][xtrack] /= n;
205 else
206 pert_15mu_low->bt[track0 + track][xtrack] = GSL_NAN;
207 }
208
209 /* Get 15 micron brightness temperature (high altitudes)... */
210 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
211 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
212 n = 0;
213 for (int i = 0; i < N15_HIGH; i++)
214 if (gsl_finite(airs_rad_gran.radiances
215 [track][xtrack][list_15mu_high[i]])) {
216 pert_15mu_high->bt[track0 + track][xtrack]
217 += BRIGHT(airs_rad_gran.radiances
218 [track][xtrack][list_15mu_high[i]],
219 airs_rad_gran.nominal_freq[list_15mu_high[i]]);
220 n++;
221 }
222 if (n > 0.9 * N15_HIGH)
223 pert_15mu_high->bt[track0 + track][xtrack] /= n;
224 else
225 pert_15mu_high->bt[track0 + track][xtrack] = GSL_NAN;
226 }
227
228 /* Increment track counter... */
229 track0 += AIRS_RAD_GEOTRACK;
230 }
231
232 /* ------------------------------------------------------------
233 Calculate perturbations and variances...
234 ------------------------------------------------------------ */
235
236 /* Convert to wave analysis struct... */
237 pert2wave(pert_4mu, &wave,
238 0, pert_4mu->ntrack - 1, 0, pert_4mu->nxtrack - 1);
239
240 /* Estimate background... */
241 background_poly(&wave, 5, 0);
242
243 /* Compute variance... */
244 variance(&wave, var_dh);
245
246 /* Copy data... */
247 for (int ix = 0; ix < wave.nx; ix++)
248 for (int iy = 0; iy < wave.ny; iy++) {
249 pert_4mu->pt[iy][ix] = wave.pt[ix][iy];
250 pert_4mu->var[iy][ix] = wave.var[ix][iy];
251 }
252
253 /* Convert to wave analysis struct... */
254 pert2wave(pert_15mu_low, &wave,
255 0, pert_15mu_low->ntrack - 1, 0, pert_15mu_low->nxtrack - 1);
256
257 /* Estimate background... */
258 background_poly(&wave, 5, 0);
259
260 /* Compute variance... */
261 variance(&wave, var_dh);
262
263 /* Copy data... */
264 for (int ix = 0; ix < wave.nx; ix++)
265 for (int iy = 0; iy < wave.ny; iy++) {
266 pert_15mu_low->pt[iy][ix] = wave.pt[ix][iy];
267 pert_15mu_low->var[iy][ix] = wave.var[ix][iy];
268 }
269
270 /* Convert to wave analysis struct... */
271 pert2wave(pert_15mu_high, &wave,
272 0, pert_15mu_high->ntrack - 1, 0, pert_15mu_high->nxtrack - 1);
273
274 /* Estimate background... */
275 background_poly(&wave, 5, 0);
276
277 /* Compute variance... */
278 variance(&wave, var_dh);
279
280 /* Copy data... */
281 for (int ix = 0; ix < wave.nx; ix++)
282 for (int iy = 0; iy < wave.ny; iy++) {
283 pert_15mu_high->pt[iy][ix] = wave.pt[ix][iy];
284 pert_15mu_high->var[iy][ix] = wave.var[ix][iy];
285 }
286
287 /* ------------------------------------------------------------
288 Write to netCDF file...
289 ------------------------------------------------------------ */
290
291 /* Create netCDF file... */
292 NC(nc_create(argv[1], NC_CLOBBER, &ncid));
293
294 /* Set dimensions... */
295 NC(nc_def_dim(ncid, "NTRACK", NC_UNLIMITED, &dimid[0]));
296 NC(nc_def_dim(ncid, "NXTRACK", AIRS_RAD_GEOXTRACK, &dimid[1]));
297
298 /* Add variables... */
299 NC(nc_def_var(ncid, "time", NC_DOUBLE, 2, dimid, &time_varid));
300 add_att(ncid, time_varid, "s", "time (seconds since 2000-01-01T00:00Z)");
301 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 2, dimid, &lon_varid));
302 add_att(ncid, lon_varid, "deg", "footprint longitude");
303 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 2, dimid, &lat_varid));
304 add_att(ncid, lat_varid, "deg", "footprint latitude");
305
306 NC(nc_def_var(ncid, "bt_8mu", NC_FLOAT, 2, dimid, &bt_8mu_varid));
307 add_att(ncid, bt_8mu_varid, "K", "brightness temperature at 8.1 micron");
308
309 NC(nc_def_var(ncid, "bt_4mu", NC_FLOAT, 2, dimid, &bt_4mu_varid));
310 add_att(ncid, bt_4mu_varid, "K", "brightness temperature" " at 4.3 micron");
311 NC(nc_def_var(ncid, "bt_4mu_pt", NC_FLOAT, 2, dimid, &bt_4mu_pt_varid));
312 add_att(ncid, bt_4mu_pt_varid, "K", "brightness temperature perturbation"
313 " at 4.3 micron");
314 NC(nc_def_var(ncid, "bt_4mu_var", NC_FLOAT, 2, dimid, &bt_4mu_var_varid));
315 add_att(ncid, bt_4mu_var_varid, "K^2", "brightness temperature variance"
316 " at 4.3 micron");
317
318 NC(nc_def_var(ncid, "bt_15mu_low", NC_FLOAT, 2, dimid, &bt_15mu_low_varid));
319 add_att(ncid, bt_15mu_low_varid, "K", "brightness temperature"
320 " at 15 micron (low altitudes)");
321 NC(nc_def_var(ncid, "bt_15mu_low_pt", NC_FLOAT, 2, dimid,
322 &bt_15mu_low_pt_varid));
323 add_att(ncid, bt_15mu_low_pt_varid, "K",
324 "brightness temperature perturbation"
325 " at 15 micron (low altitudes)");
326 NC(nc_def_var
327 (ncid, "bt_15mu_low_var", NC_FLOAT, 2, dimid, &bt_15mu_low_var_varid));
328 add_att(ncid, bt_15mu_low_var_varid, "K^2",
329 "brightness temperature variance" " at 15 micron (low altitudes)");
330
331 NC(nc_def_var(ncid, "bt_15mu_high", NC_FLOAT, 2, dimid,
332 &bt_15mu_high_varid));
333 add_att(ncid, bt_15mu_high_varid, "K", "brightness temperature"
334 " at 15 micron (high altitudes)");
335 NC(nc_def_var(ncid, "bt_15mu_high_pt", NC_FLOAT, 2, dimid,
336 &bt_15mu_high_pt_varid));
337 add_att(ncid, bt_15mu_high_pt_varid, "K",
338 "brightness temperature perturbation"
339 " at 15 micron (high altitudes)");
340 NC(nc_def_var
341 (ncid, "bt_15mu_high_var", NC_FLOAT, 2, dimid, &bt_15mu_high_var_varid));
342 add_att(ncid, bt_15mu_high_var_varid, "K^2",
343 "brightness temperature variance" " at 15 micron (high altitudes)");
344
345 /* Leave define mode... */
346 NC(nc_enddef(ncid));
347
348 /* Loop over tracks... */
349 for (int track = 0; track < pert_4mu->ntrack; track++) {
350
351 /* Set array sizes... */
352 start[0] = (size_t) track;
353 start[1] = 0;
354 count[0] = 1;
355 count[1] = (size_t) pert_4mu->nxtrack;
356
357 /* Write data... */
358 NC(nc_put_vara_double(ncid, time_varid, start, count,
359 pert_4mu->time[track]));
360 NC(nc_put_vara_double(ncid, lon_varid, start, count,
361 pert_4mu->lon[track]));
362 NC(nc_put_vara_double(ncid, lat_varid, start, count,
363 pert_4mu->lat[track]));
364
365 NC(nc_put_vara_double(ncid, bt_8mu_varid, start, count,
366 pert_4mu->dc[track]));
367
368 NC(nc_put_vara_double(ncid, bt_4mu_varid, start, count,
369 pert_4mu->bt[track]));
370 NC(nc_put_vara_double(ncid, bt_4mu_pt_varid, start, count,
371 pert_4mu->pt[track]));
372 NC(nc_put_vara_double(ncid, bt_4mu_var_varid, start, count,
373 pert_4mu->var[track]));
374
375 NC(nc_put_vara_double(ncid, bt_15mu_low_varid, start, count,
376 pert_15mu_low->bt[track]));
377 NC(nc_put_vara_double(ncid, bt_15mu_low_pt_varid, start, count,
378 pert_15mu_low->pt[track]));
379 NC(nc_put_vara_double(ncid, bt_15mu_low_var_varid, start, count,
380 pert_15mu_low->var[track]));
381
382 NC(nc_put_vara_double(ncid, bt_15mu_high_varid, start, count,
383 pert_15mu_high->bt[track]));
384 NC(nc_put_vara_double(ncid, bt_15mu_high_pt_varid, start, count,
385 pert_15mu_high->pt[track]));
386 NC(nc_put_vara_double(ncid, bt_15mu_high_var_varid, start, count,
387 pert_15mu_high->var[track]));
388 }
389
390 /* Close file... */
391 NC(nc_close(ncid));
392
393 /* Free... */
394 free(pert_4mu);
395 free(pert_15mu_low);
396 free(pert_15mu_high);
397
398 return EXIT_SUCCESS;
399}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: diff_apr.c:35
BRIGHT
#define BRIGHT(rad, nu)
Compute brightness temperature.
Definition: jurassic.h:126
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
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
libairs.h
AIRS Code Collection library declarations.
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libairs.h:126
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libairs.h:123
main
int main(int argc, char *argv[])
Definition: perturbation.c:45
N15_LOW
#define N15_LOW
Definition: perturbation.c:36
N4
#define N4
Definition: perturbation.c:33
N15_HIGH
#define N15_HIGH
Definition: perturbation.c:39
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
Generated on Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4