AIRS Code Collection
volcano.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 Functions...
30 ------------------------------------------------------------ */
31
32/* Estimate noise. */
33double get_noise(
34 double bt,
35 double dt250,
36 double nu);
37
38/* ------------------------------------------------------------
39 Main...
40 ------------------------------------------------------------ */
41
42int main(
43 int argc,
44 char *argv[]) {
45
46 FILE *out;
47
48 static airs_rad_gran_t airs_rad_gran;
49
50 const double ci_nedt = 0.0783, ai_low_bt1_nedt = 0.3698, ai_low_bt2_nedt =
51 0.1177, ai_high_bt1_nedt = 0.0766, ai_high_bt2_nedt =
52 0.3706, ai_old_bt1_nedt = 0.3155, ai_old_bt2_nedt =
53 0.1177, si_high_bt1_nedt = 0.1025, si_high_bt2_nedt =
54 0.1373, si_low_bt1_nedt = 0.0799, si_low_bt2_nedt =
55 0.0909, si_old_bt1_nedt = 0.1064, si_old_bt2_nedt =
56 0.0909, si_oper_bt1_nedt = 0.0884, si_oper_bt2_nedt = 0.1159;
57
58 const int ai_low_nu1 = 641, ai_low_nu2 = 901, ai_high_nu1 =
59 1295, ai_high_nu2 = 1162, ai_old_nu1 = 559, ai_old_nu2 = 901, ci_nu =
60 1290, si_low_nu1 = 1601, si_low_nu2 = 1526, si_high_nu1 =
61 1602, si_high_nu2 = 1551, si_old_nu1 = 1591, si_old_nu2 =
62 1526, si_oper_nu1 = 1636, si_oper_nu2 = 1507;
63
64 /* Check arguments... */
65 if (argc < 3)
66 ERRMSG("Give parameters: <out.tab> <l1b_file1> [<l1b_file2> ...]");
67
68 /* Create file... */
69 printf("Write volcanic emission data: %s\n", argv[1]);
70 if (!(out = fopen(argv[1], "w")))
71 ERRMSG("Cannot create file!");
72
73 /* Loop over HDF files... */
74 for (int iarg = 2; iarg < argc; iarg++) {
75
76 /* Read AIRS data... */
77 printf("Read AIRS Level-1B data file: %s\n", argv[iarg]);
78 airs_rad_rdr(argv[iarg], &airs_rad_gran);
79
80 /* Write header... */
81 if (iarg == 2) {
82 fprintf(out,
83 "# $1 = time [s]\n"
84 "# $2 = footprint longitude [deg]\n"
85 "# $3 = footprint latitude [deg]\n"
86 "# $4 = satellite altitude [km]\n"
87 "# $5 = satellite longitude [deg]\n"
88 "# $6 = satellite latitude [deg]\n");
89 fprintf(out,
90 "# $7 = cloud index, BT(%.2f/cm) [K]\n"
91 "# $8 = cloud index error [K]\n"
92 "# $9 = ash index (low wavenumbers),"
93 " BT(%.2f/cm) - BT(%.2f/cm) [K]\n"
94 "# $10 = ash index (low wavenumbers) error [K]\n"
95 "# $11 = ash index (high wavenumbers),"
96 " BT(%.2f/cm) - BT(%.2f/cm) [K]\n"
97 "# $12 = ash index (high wavenumbers) error [K]\n"
98 "# $13 = ash index (Hoffmann et al., 2014),"
99 " BT(%.2f/cm) - BT(%.2f/cm) [K]\n"
100 "# $14 = ash index (Hoffmann et al., 2014) error [K]\n",
101 airs_rad_gran.nominal_freq[ci_nu],
102 airs_rad_gran.nominal_freq[ai_low_nu1],
103 airs_rad_gran.nominal_freq[ai_low_nu2],
104 airs_rad_gran.nominal_freq[ai_high_nu1],
105 airs_rad_gran.nominal_freq[ai_high_nu2],
106 airs_rad_gran.nominal_freq[ai_old_nu1],
107 airs_rad_gran.nominal_freq[ai_old_nu2]);
108 fprintf(out,
109 "# $15 = SO2 index (low concentrations),"
110 " BT(%.2f/cm) - BT(%.2f/cm) [K]\n"
111 "# $16 = SO2 index (low concentrations) error [K]\n"
112 "# $17 = SO2 index (high concentrations),"
113 " BT(%.2f/cm) - BT(%.2f/cm) [K]\n"
114 "# $18 = SO2 index (high concentrations) error [K]\n"
115 "# $19 = SO2 index (operational),"
116 " BT(%.2f/cm) - BT(%.2f/cm) [K]\n"
117 "# $20 = SO2 index (operational) error [K]\n"
118 "# $21 = SO2 index (Hoffmann et al., 2014),"
119 " BT(%.2f/cm) - BT(%.2f/cm) [K]\n"
120 "# $22 = SO2 index (Hoffmann et al., 2014) error [K]\n",
121 airs_rad_gran.nominal_freq[si_low_nu1],
122 airs_rad_gran.nominal_freq[si_low_nu2],
123 airs_rad_gran.nominal_freq[si_high_nu1],
124 airs_rad_gran.nominal_freq[si_high_nu2],
125 airs_rad_gran.nominal_freq[si_oper_nu1],
126 airs_rad_gran.nominal_freq[si_oper_nu2],
127 airs_rad_gran.nominal_freq[si_old_nu1],
128 airs_rad_gran.nominal_freq[si_old_nu2]);
129 }
130
131 /* Flag bad observations... */
132 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
133 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++)
134 for (int ichan = 0; ichan < AIRS_RAD_CHANNEL; ichan++)
135 if ((airs_rad_gran.state[track][xtrack] != 0)
136 || (airs_rad_gran.ExcludedChans[ichan] > 2)
137 || (airs_rad_gran.CalChanSummary[ichan] & 8)
138 || (airs_rad_gran.CalChanSummary[ichan] & (32 + 64))
139 || (airs_rad_gran.CalFlag[track][ichan] & 16))
140 airs_rad_gran.radiances[track][xtrack][ichan] = GSL_NAN;
141
142 /* Loop over scans... */
143 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++) {
144
145 /* Write output... */
146 fprintf(out, "\n");
147
148 /* Loop over footprints... */
149 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
150
151 /* cloud index... */
152 double ci =
153 BRIGHT(airs_rad_gran.radiances[track][xtrack][ci_nu] * 0.001,
154 airs_rad_gran.nominal_freq[ci_nu]);
155 double ci_err =
156 get_noise(ci, ci_nedt, airs_rad_gran.nominal_freq[ci_nu]);
157
158 /* ash index (low wavenumbers)... */
159 double ai_low_bt1 =
160 BRIGHT(airs_rad_gran.radiances[track][xtrack][ai_low_nu1] *
161 0.001, airs_rad_gran.nominal_freq[ai_low_nu1]);
162 double ai_low_bt2 =
163 BRIGHT(airs_rad_gran.radiances[track][xtrack][ai_low_nu2] *
164 0.001, airs_rad_gran.nominal_freq[ai_low_nu2]);
165 double ai_low = ai_low_bt1 - ai_low_bt2;
166 double ai_low_err =
167 sqrt(gsl_pow_2(get_noise(ai_low_bt1, ai_low_bt1_nedt,
168 airs_rad_gran.nominal_freq[ai_low_nu1]))
169 + gsl_pow_2(get_noise(ai_low_bt2, ai_low_bt2_nedt,
170 airs_rad_gran.nominal_freq
171 [ai_low_nu2])));
172
173 /* ash index (high wavenumbers)... */
174 double ai_high_bt1 =
175 BRIGHT(airs_rad_gran.radiances[track][xtrack][ai_high_nu1] *
176 0.001, airs_rad_gran.nominal_freq[ai_high_nu1]);
177 double ai_high_bt2 =
178 BRIGHT(airs_rad_gran.radiances[track][xtrack][ai_high_nu2] *
179 0.001, airs_rad_gran.nominal_freq[ai_high_nu2]);
180 double ai_high = ai_high_bt1 - ai_high_bt2;
181 double ai_high_err =
182 sqrt(gsl_pow_2(get_noise(ai_high_bt1, ai_high_bt1_nedt,
183 airs_rad_gran.nominal_freq[ai_high_nu1]))
184 + gsl_pow_2(get_noise(ai_high_bt2, ai_high_bt2_nedt,
185 airs_rad_gran.nominal_freq
186 [ai_high_nu2])));
187
188 /* ash index (old)... */
189 double ai_old_bt1 =
190 BRIGHT(airs_rad_gran.radiances[track][xtrack][ai_old_nu1] *
191 0.001, airs_rad_gran.nominal_freq[ai_old_nu1]);
192 double ai_old_bt2 =
193 BRIGHT(airs_rad_gran.radiances[track][xtrack][ai_old_nu2] *
194 0.001, airs_rad_gran.nominal_freq[ai_old_nu2]);
195 double ai_old = ai_old_bt1 - ai_old_bt2;
196 double ai_old_err =
197 sqrt(gsl_pow_2(get_noise(ai_old_bt1, ai_old_bt1_nedt,
198 airs_rad_gran.nominal_freq[ai_old_nu1]))
199 + gsl_pow_2(get_noise(ai_old_bt2, ai_old_bt2_nedt,
200 airs_rad_gran.nominal_freq
201 [ai_old_nu2])));
202
203 /* SO2 index (low concentrations)... */
204 double si_low_bt1 =
205 BRIGHT(airs_rad_gran.radiances[track][xtrack][si_low_nu1] *
206 0.001, airs_rad_gran.nominal_freq[si_low_nu1]);
207 double si_low_bt2 =
208 BRIGHT(airs_rad_gran.radiances[track][xtrack][si_low_nu2] *
209 0.001, airs_rad_gran.nominal_freq[si_low_nu2]);
210 double si_low = si_low_bt1 - si_low_bt2;
211 double si_low_err =
212 sqrt(gsl_pow_2(get_noise(si_low_bt1, si_low_bt1_nedt,
213 airs_rad_gran.nominal_freq[si_low_nu1]))
214 + gsl_pow_2(get_noise(si_low_bt2, si_low_bt2_nedt,
215 airs_rad_gran.nominal_freq
216 [si_low_nu2])));
217
218 /* SO2 index (high concentrations)... */
219 double si_high_bt1 =
220 BRIGHT(airs_rad_gran.radiances[track][xtrack][si_high_nu1] *
221 0.001, airs_rad_gran.nominal_freq[si_high_nu1]);
222 double si_high_bt2 =
223 BRIGHT(airs_rad_gran.radiances[track][xtrack][si_high_nu2] *
224 0.001, airs_rad_gran.nominal_freq[si_high_nu2]);
225 double si_high = si_high_bt1 - si_high_bt2;
226 double si_high_err =
227 sqrt(gsl_pow_2(get_noise(si_high_bt1, si_high_bt1_nedt,
228 airs_rad_gran.nominal_freq[si_high_nu1]))
229 + gsl_pow_2(get_noise(si_high_bt2, si_high_bt2_nedt,
230 airs_rad_gran.nominal_freq
231 [si_high_nu2])));
232
233 /* SO2 index (operational)... */
234 double si_oper_bt1 =
235 BRIGHT(airs_rad_gran.radiances[track][xtrack][si_oper_nu1] *
236 0.001, airs_rad_gran.nominal_freq[si_oper_nu1]);
237 double si_oper_bt2 =
238 BRIGHT(airs_rad_gran.radiances[track][xtrack][si_oper_nu2] *
239 0.001, airs_rad_gran.nominal_freq[si_oper_nu2]);
240 double si_oper = si_oper_bt1 - si_oper_bt2;
241 double si_oper_err =
242 sqrt(gsl_pow_2(get_noise(si_oper_bt1, si_oper_bt1_nedt,
243 airs_rad_gran.nominal_freq[si_oper_nu1]))
244 + gsl_pow_2(get_noise(si_oper_bt2, si_oper_bt2_nedt,
245 airs_rad_gran.nominal_freq
246 [si_oper_nu2])));
247
248 /* SO2 index (old)... */
249 double si_old_bt1 =
250 BRIGHT(airs_rad_gran.radiances[track][xtrack][si_old_nu1] *
251 0.001, airs_rad_gran.nominal_freq[si_old_nu1]);
252 double si_old_bt2 =
253 BRIGHT(airs_rad_gran.radiances[track][xtrack][si_old_nu2] *
254 0.001, airs_rad_gran.nominal_freq[si_old_nu2]);
255 double si_old = si_old_bt1 - si_old_bt2;
256 double si_old_err =
257 sqrt(gsl_pow_2(get_noise(si_old_bt1, si_old_bt1_nedt,
258 airs_rad_gran.nominal_freq[si_old_nu1]))
259 + gsl_pow_2(get_noise(si_old_bt2, si_old_bt2_nedt,
260 airs_rad_gran.nominal_freq
261 [si_old_nu2])));
262
263 /* Write output... */
264 fprintf(out,
265 "%.2f %.4f %.4f %.3f %.4f %.4f %.2f %.2f %.2f %.2f %.2f %.2f "
266 "%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",
267 airs_rad_gran.Time[track][xtrack] - 220838400,
268 airs_rad_gran.Longitude[track][xtrack],
269 airs_rad_gran.Latitude[track][xtrack],
270 airs_rad_gran.satheight[track],
271 airs_rad_gran.sat_lon[track],
272 airs_rad_gran.sat_lat[track],
273 ci, ci_err, ai_low, ai_low_err, ai_high, ai_high_err, ai_old,
274 ai_old_err, si_low, si_low_err, si_high, si_high_err, si_oper,
275 si_oper_err, si_old, si_old_err);
276 }
277 }
278 }
279
280 /* Close file... */
281 fclose(out);
282
283 return EXIT_SUCCESS;
284}
285
286/************************************************************************/
287
288double get_noise(
289 double bt,
290 double dt250,
291 double nu) {
292
293 double nesr;
294
295 nesr = PLANCK(250.0 + dt250, nu) - PLANCK(250.0, nu);
296
297 return BRIGHT(PLANCK(bt, nu) + nesr, nu) - bt;
298}
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
PLANCK
#define PLANCK(T, nu)
Compute Planck function.
Definition: jurassic.h:183
libairs.h
AIRS Code Collection library declarations.
main
int main(int argc, char *argv[])
Definition: volcano.c:42
get_noise
double get_noise(double bt, double dt250, double nu)
Definition: volcano.c:288
Generated on Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4