AIRS Code Collection
nlte.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 <mpi.h>
27#include <omp.h>
28#include <netcdf.h>
29#include "jurassic.h"
30
31/* ------------------------------------------------------------
32 Macros...
33 ------------------------------------------------------------ */
34
36#define NC(cmd) { \
37 int nc_result=(cmd); \
38 if(nc_result!=NC_NOERR) \
39 ERRMSG("%s", nc_strerror(nc_result)); \
40}
41
42/* ------------------------------------------------------------
43 Dimensions...
44 ------------------------------------------------------------ */
45
47#define L1_NCHAN 34
48
50#define L1_NTRACK 135
51
53#define L1_NXTRACK 90
54
56#define L2_NLAY 27
57
59#define L2_NTRACK 45
60
62#define L2_NXTRACK 30
63
64/* ------------------------------------------------------------
65 Structs...
66 ------------------------------------------------------------ */
67
69typedef struct {
70
72 int ncid;
73
75 int np;
76
78 double l1_time[L1_NTRACK][L1_NXTRACK];
79
81 double l1_lon[L1_NTRACK][L1_NXTRACK];
82
84 double l1_lat[L1_NTRACK][L1_NXTRACK];
85
87 double l1_sat_z[L1_NTRACK];
88
90 double l1_sat_lon[L1_NTRACK];
91
93 double l1_sat_lat[L1_NTRACK];
94
96 double l1_nu[L1_NCHAN];
97
99 float l1_rad[L1_NTRACK][L1_NXTRACK][L1_NCHAN];
100
102 double l2_z[L2_NTRACK][L2_NXTRACK][L2_NLAY];
103
105 double l2_p[L2_NLAY];
106
108 double l2_t[L2_NTRACK][L2_NXTRACK][L2_NLAY];
109
111 float ret_z[NP];
112
114 float ret_p[L1_NTRACK * L1_NXTRACK];
115
117 float ret_t[L1_NTRACK * L1_NXTRACK * NP];
118
119} ncd_t;
120
121/* ------------------------------------------------------------
122 Functions...
123 ------------------------------------------------------------ */
124
126void add_var(
127 int ncid,
128 const char *varname,
129 const char *unit,
130 const char *longname,
131 int type,
132 int dimid[],
133 int *varid,
134 int ndims);
135
137void fill_gaps(
138 double x[L2_NTRACK][L2_NXTRACK][L2_NLAY],
139 double cx,
140 double cy);
141
143void init_l2(
144 ncd_t * ncd,
145 int track,
146 int xtrack,
147 ctl_t * ctl,
148 atm_t * atm);
149
151void read_nc(
152 char *filename,
153 ncd_t * ncd);
154
155/* ------------------------------------------------------------
156 Main...
157 ------------------------------------------------------------ */
158
159int main(
160 int argc,
161 char *argv[]) {
162
163 static ctl_t ctl;
164 static atm_t atm_apr, atm_clim, atm_i;
165 static obs_t obs_i, obs_meas;
166 static ncd_t ncd;
167 static ret_t ret;
168
169 static FILE *in, *out;
170
171 static char filename[LEN], filename2[2 * LEN];
172
173 static double chisq[L1_NTRACK][L1_NXTRACK], ni[L1_NTRACK][L1_NXTRACK],
174 z[NP];
175
176 static int channel[ND], n, ntask = -1, rank, size;
177
178 /* ------------------------------------------------------------
179 Init...
180 ------------------------------------------------------------ */
181
182 /* MPI... */
183 MPI_Init(&argc, &argv);
184 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
185 MPI_Comm_size(MPI_COMM_WORLD, &size);
186
187 /* Measure CPU time... */
188 TIMER("total", 1);
189
190 /* Check arguments... */
191 if (argc < 3)
192 ERRMSG("Give parameters: <ctl> <filelist>");
193
194 /* Read control parameters... */
195 read_ctl(argc, argv, &ctl);
196 read_ret(argc, argv, &ctl, &ret);
197
198 /* Initialize look-up tables... */
199 tbl_t *tbl = read_tbl(&ctl);
200
201 /* Read retrieval grid... */
202 const int nz = (int) scan_ctl(argc, argv, "NZ", -1, "", NULL);
203 if (nz > NP)
204 ERRMSG("Too many altitudes!");
205 for (int iz = 0; iz < nz; iz++)
206 z[iz] = scan_ctl(argc, argv, "Z", iz, "", NULL);
207
208 /* Read track range... */
209 const int track0 = (int) scan_ctl(argc, argv, "TRACK_MIN", -1, "0", NULL);
210 const int track1 = (int) scan_ctl(argc, argv, "TRACK_MAX", -1, "134", NULL);
211
212 /* Read xtrack range... */
213 const int xtrack0 = (int) scan_ctl(argc, argv, "XTRACK_MIN", -1, "0", NULL);
214 const int xtrack1 =
215 (int) scan_ctl(argc, argv, "XTRACK_MAX", -1, "89", NULL);
216
217 /* Background smoothing... */
218 const double sx = scan_ctl(argc, argv, "SX", -1, "8", NULL);
219 const double sy = scan_ctl(argc, argv, "SY", -1, "2", NULL);
220
221 /* ------------------------------------------------------------
222 Distribute granules...
223 ------------------------------------------------------------ */
224
225 /* Open filelist... */
226 printf("Read filelist: %s\n", argv[2]);
227 if (!(in = fopen(argv[2], "r")))
228 ERRMSG("Cannot open filelist!");
229
230 /* Loop over netCDF files... */
231 while (fscanf(in, "%s", filename) != EOF) {
232
233 /* Distribute files with MPI... */
234 if ((++ntask) % size != rank)
235 continue;
236
237 /* Write info... */
238 printf("Retrieve file %s on rank %d of %d (with %d threads)...\n",
239 filename, rank + 1, size, omp_get_max_threads());
240
241 /* ------------------------------------------------------------
242 Initialize retrieval...
243 ------------------------------------------------------------ */
244
245 /* Read netCDF file... */
246 read_nc(filename, &ncd);
247
248 /* Identify radiance channels... */
249 for (int id = 0; id < ctl.nd; id++) {
250 channel[id] = -999;
251 for (int i = 0; i < L1_NCHAN; i++)
252 if (fabs(ctl.nu[id] - ncd.l1_nu[i]) < 0.1)
253 channel[id] = i;
254 if (channel[id] < 0)
255 ERRMSG("Cannot identify radiance channel!");
256 }
257
258 /* Fill data gaps... */
259 fill_gaps(ncd.l2_t, sx, sy);
260 fill_gaps(ncd.l2_z, sx, sy);
261
262 /* Set climatological data for center of granule... */
263 atm_clim.np = nz;
264 for (int iz = 0; iz < nz; iz++)
265 atm_clim.z[iz] = z[iz];
266 climatology(&ctl, &atm_clim);
267
268 /* ------------------------------------------------------------
269 Retrieval...
270 ------------------------------------------------------------ */
271
272 /* Loop over swaths... */
273 for (int track = track0; track <= track1; track++) {
274
275 /* Measure CPU time... */
276 TIMER("retrieval", 1);
277
278 /* Loop over scan... */
279 for (int xtrack = xtrack0; xtrack <= xtrack1; xtrack++) {
280
281 /* Store observation data... */
282 obs_meas.nr = 1;
283 obs_meas.time[0] = ncd.l1_time[track][xtrack];
284 obs_meas.obsz[0] = ncd.l1_sat_z[track];
285 obs_meas.obslon[0] = ncd.l1_sat_lon[track];
286 obs_meas.obslat[0] = ncd.l1_sat_lat[track];
287 obs_meas.vplon[0] = ncd.l1_lon[track][xtrack];
288 obs_meas.vplat[0] = ncd.l1_lat[track][xtrack];
289 for (int id = 0; id < ctl.nd; id++)
290 obs_meas.rad[id][0] = ncd.l1_rad[track][xtrack][channel[id]];
291
292 /* Flag out 4 micron channels... */
293 for (int id = 0; id < ctl.nd; id++)
294 if (ctl.nu[id] >= 2000)
295 obs_meas.rad[id][0] = GSL_NAN;
296
297 /* Prepare atmospheric data... */
298 copy_atm(&ctl, &atm_apr, &atm_clim, 0);
299 for (int ip = 0; ip < atm_apr.np; ip++) {
300 atm_apr.time[ip] = obs_meas.time[0];
301 atm_apr.lon[ip] = obs_meas.vplon[0];
302 atm_apr.lat[ip] = obs_meas.vplat[0];
303 }
304
305 /* Merge Level-2 data... */
306 init_l2(&ncd, track, xtrack, &ctl, &atm_apr);
307
308 /* Retrieval... */
309 optimal_estimation(&ret, &ctl, tbl, &obs_meas, &obs_i,
310 &atm_apr, &atm_i, &chisq[track][xtrack]);
311
312 /* Run forward model including 4 micron channels... */
313 for (int id = 0; id < ctl.nd; id++)
314 obs_meas.rad[id][0] = obs_i.rad[id][0]
315 = ncd.l1_rad[track][xtrack][channel[id]];
316 formod(&ctl, tbl, &atm_i, &obs_i);
317
318 /* Calculate non-LTE index... */
319 n = 0;
320 ni[track][xtrack] = 0;
321 for (int id = 0; id < ctl.nd; id++)
322 if (ctl.nu[id] >= 2000 && gsl_finite(obs_meas.rad[id][0])) {
323 ni[track][xtrack] +=
324 (BRIGHT(obs_meas.rad[id][0], ncd.l1_nu[channel[id]])
325 - BRIGHT(obs_i.rad[id][0], ncd.l1_nu[channel[id]]));
326 n++;
327 }
328 ni[track][xtrack] /= n;
329 }
330
331 /* Measure CPU time... */
332 TIMER("retrieval", 3);
333 }
334
335 /* ------------------------------------------------------------
336 Write output...
337 ------------------------------------------------------------ */
338
339 /* Set filename... */
340 sprintf(filename2, "%s.nlte.tab", filename);
341
342 /* Create output file... */
343 printf("Write non-LTE data: %s\n", filename2);
344 if (!(out = fopen(filename2, "w")))
345 ERRMSG("Cannot create file!");
346
347 /* Write header... */
348 fprintf(out,
349 "# $1 = time (seconds since 2000-01-01, 00:00 UTC)\n"
350 "# $2 = longitude [deg]\n"
351 "# $3 = latitude [deg]\n"
352 "# $4 = solar zenith angle [deg]\n"
353 "# $5 = non-LTE index [K]\n" "# $6 = chi^2 of retrieval fit\n");
354
355 /* Write data... */
356 for (int track = track0; track <= track1; track++) {
357 fprintf(out, "\n");
358 for (int xtrack = xtrack0; xtrack <= xtrack1; xtrack++)
359 fprintf(out, "%.2f %g %g %g %g %g\n",ncd.l1_time[track][xtrack],
360 ncd.l1_lon[track][xtrack], ncd.l1_lat[track][xtrack],
361 RAD2DEG(acos(cos_sza(ncd.l1_time[track][xtrack],
362 ncd.l1_lon[track][xtrack],
363 ncd.l1_lat[track][xtrack]))),
364 ni[track][xtrack], chisq[track][xtrack]);
365 }
366
367 /* Close output file... */
368 fclose(out);
369
370 /* Write info... */
371 printf("Retrieval finished on rank %d of %d!\n", rank, size);
372 }
373
374 /* Close file list... */
375 fclose(in);
376
377 /* Free... */
378 free(tbl);
379
380 /* Measure CPU time... */
381 TIMER("total", 3);
382
383 /* Report memory usage... */
384 printf("MEMORY_ATM = %g MByte\n", 4. * sizeof(atm_t) / 1024. / 1024.);
385 printf("MEMORY_CTL = %g MByte\n", 1. * sizeof(ctl_t) / 1024. / 1024.);
386 printf("MEMORY_NCD = %g MByte\n", 1. * sizeof(ncd_t) / 1024. / 1024.);
387 printf("MEMORY_OBS = %g MByte\n", 3. * sizeof(atm_t) / 1024. / 1024.);
388 printf("MEMORY_RET = %g MByte\n", 1. * sizeof(ret_t) / 1024. / 1024.);
389 printf("MEMORY_TBL = %g MByte\n", 1. * sizeof(tbl_t) / 1024. / 1024.);
390
391 /* Report problem size... */
392 printf("SIZE_TASKS = %d\n", size);
393 printf("SIZE_THREADS = %d\n", omp_get_max_threads());
394
395 /* MPI... */
396 MPI_Finalize();
397
398 return EXIT_SUCCESS;
399}
400
401/*****************************************************************************/
402
403void add_var(
404 int ncid,
405 const char *varname,
406 const char *unit,
407 const char *longname,
408 int type,
409 int dimid[],
410 int *varid,
411 int ndims) {
412
413 /* Check if variable exists... */
414 if (nc_inq_varid(ncid, varname, varid) != NC_NOERR) {
415
416 /* Define variable... */
417 NC(nc_def_var(ncid, varname, type, ndims, dimid, varid));
418
419 /* Set long name... */
420 NC(nc_put_att_text
421 (ncid, *varid, "long_name", strlen(longname), longname));
422
423 /* Set units... */
424 NC(nc_put_att_text(ncid, *varid, "units", strlen(unit), unit));
425 }
426}
427
428/************************************************************************/
429
430void fill_gaps(
431 double x[L2_NTRACK][L2_NXTRACK][L2_NLAY],
432 double cx,
433 double cy) {
434
435 double help[L2_NTRACK][L2_NXTRACK];
436
437 /* Loop over layers... */
438 for (int lay = 0; lay < L2_NLAY; lay++) {
439
440 /* Loop over grid points... */
441 for (int track = 0; track < L2_NTRACK; track++)
442 for (int xtrack = 0; xtrack < L2_NXTRACK; xtrack++) {
443
444 /* Init... */
445 help[track][xtrack] = 0;
446 double wsum = 0;
447
448 /* Averrage data points... */
449 for (int track2 = 0; track2 < L2_NTRACK; track2++)
450 for (int xtrack2 = 0; xtrack2 < L2_NXTRACK; xtrack2++)
451 if (gsl_finite(x[track2][xtrack2][lay])
452 && x[track2][xtrack2][lay] > 0) {
453 const double w = exp(-gsl_pow_2((xtrack - xtrack2) / cx)
454 - gsl_pow_2((track - track2) / cy));
455 help[track][xtrack] += w * x[track2][xtrack2][lay];
456 wsum += w;
457 }
458
459 /* Normalize... */
460 if (wsum > 0)
461 help[track][xtrack] /= wsum;
462 else
463 help[track][xtrack] = GSL_NAN;
464 }
465
466 /* Copy grid points... */
467 for (int track = 0; track < L2_NTRACK; track++)
468 for (int xtrack = 0; xtrack < L2_NXTRACK; xtrack++)
469 x[track][xtrack][lay] = help[track][xtrack];
470 }
471}
472
473/************************************************************************/
474
475void init_l2(
476 ncd_t *ncd,
477 int track,
478 int xtrack,
479 ctl_t *ctl,
480 atm_t *atm) {
481
482 static atm_t atm_airs;
483
484 double k[NW], p, q[NG], t, zmax = 0, zmin = 1000;
485
486 /* Reset track- and xtrack-index to match Level-2 data... */
487 track /= 3;
488 xtrack /= 3;
489
490 /* Store AIRS data in atmospheric data struct... */
491 atm_airs.np = 0;
492 for (int lay = 0; lay < L2_NLAY; lay++)
493 if (gsl_finite(ncd->l2_z[track][xtrack][lay])) {
494 atm_airs.z[atm_airs.np] = ncd->l2_z[track][xtrack][lay];
495 atm_airs.p[atm_airs.np] = ncd->l2_p[lay];
496 atm_airs.t[atm_airs.np] = ncd->l2_t[track][xtrack][lay];
497 if ((++atm_airs.np) > NP)
498 ERRMSG("Too many layers!");
499 }
500
501 /* Check number of levels... */
502 if (atm_airs.np <= 0)
503 return;
504
505 /* Get height range of AIRS data... */
506 for (int ip = 0; ip < atm_airs.np; ip++) {
507 zmax = GSL_MAX(zmax, atm_airs.z[ip]);
508 zmin = GSL_MIN(zmin, atm_airs.z[ip]);
509 }
510
511 /* Merge AIRS data... */
512 for (int ip = 0; ip < atm->np; ip++) {
513
514 /* Interpolate AIRS data... */
515 intpol_atm(ctl, &atm_airs, atm->z[ip], &p, &t, q, k);
516
517 /* Weighting factor... */
518 double w = 1;
519 if (atm->z[ip] > zmax)
520 w = GSL_MAX(1 - (atm->z[ip] - zmax) / 50, 0);
521 if (atm->z[ip] < zmin)
522 w = GSL_MAX(1 - (zmin - atm->z[ip]) / 50, 0);
523
524 /* Merge... */
525 atm->t[ip] = w * t + (1 - w) * atm->t[ip];
526 atm->p[ip] = w * p + (1 - w) * atm->p[ip];
527 }
528}
529
530/*****************************************************************************/
531
532void read_nc(
533 char *filename,
534 ncd_t *ncd) {
535
536 int varid;
537
538 /* Open netCDF file... */
539 printf("Read netCDF file: %s\n", filename);
540 NC(nc_open(filename, NC_WRITE, &ncd->ncid));
541
542 /* Read Level-1 data... */
543 NC(nc_inq_varid(ncd->ncid, "l1_time", &varid));
544 NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_time[0]));
545 NC(nc_inq_varid(ncd->ncid, "l1_lon", &varid));
546 NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_lon[0]));
547 NC(nc_inq_varid(ncd->ncid, "l1_lat", &varid));
548 NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_lat[0]));
549 NC(nc_inq_varid(ncd->ncid, "l1_sat_z", &varid));
550 NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_sat_z));
551 NC(nc_inq_varid(ncd->ncid, "l1_sat_lon", &varid));
552 NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_sat_lon));
553 NC(nc_inq_varid(ncd->ncid, "l1_sat_lat", &varid));
554 NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_sat_lat));
555 NC(nc_inq_varid(ncd->ncid, "l1_nu", &varid));
556 NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_nu));
557 NC(nc_inq_varid(ncd->ncid, "l1_rad", &varid));
558 NC(nc_get_var_float(ncd->ncid, varid, ncd->l1_rad[0][0]));
559
560 /* Read Level-2 data... */
561 NC(nc_inq_varid(ncd->ncid, "l2_z", &varid));
562 NC(nc_get_var_double(ncd->ncid, varid, ncd->l2_z[0][0]));
563 NC(nc_inq_varid(ncd->ncid, "l2_press", &varid));
564 NC(nc_get_var_double(ncd->ncid, varid, ncd->l2_p));
565 NC(nc_inq_varid(ncd->ncid, "l2_temp", &varid));
566 NC(nc_get_var_double(ncd->ncid, varid, ncd->l2_t[0][0]));
567}
main
int main(int argc, char *argv[])
Definition: nlte.c:159
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: nlte.c:36
fill_gaps
void fill_gaps(double x[L2_NTRACK][L2_NXTRACK][L2_NLAY], double cx, double cy)
Fill data gaps in L2 data.
Definition: nlte.c:430
L1_NXTRACK
#define L1_NXTRACK
Across-track size of AIRS radiance granule (don't change).
Definition: nlte.c:53
L1_NTRACK
#define L1_NTRACK
Along-track size of AIRS radiance granule (don't change).
Definition: nlte.c:50
L2_NXTRACK
#define L2_NXTRACK
Across-track size of AIRS retrieval granule (don't change).
Definition: nlte.c:62
read_nc
void read_nc(char *filename, ncd_t *ncd)
Read netCDF file.
Definition: nlte.c:532
add_var
void add_var(int ncid, const char *varname, const char *unit, const char *longname, int type, int dimid[], int *varid, int ndims)
Create variable in netCDF file.
Definition: nlte.c:403
L2_NLAY
#define L2_NLAY
Number of AIRS pressure layers (don't change).
Definition: nlte.c:56
L1_NCHAN
#define L1_NCHAN
Number of AIRS radiance channels (don't change).
Definition: nlte.c:47
L2_NTRACK
#define L2_NTRACK
Along-track size of AIRS retrieval granule (don't change).
Definition: nlte.c:59
init_l2
void init_l2(ncd_t *ncd, int track, int xtrack, ctl_t *ctl, atm_t *atm)
Initialize with AIRS Level-2 data.
Definition: nlte.c:475
ncd_t
Buffer for netCDF data.
Definition: diff_apr.c:68
ncd_t::l1_lat
double l1_lat[L1_NTRACK][L1_NXTRACK]
Footprint latitude [deg].
Definition: diff_apr.c:83
ncd_t::l1_lon
double l1_lon[L1_NTRACK][L1_NXTRACK]
Footprint longitude [deg].
Definition: diff_apr.c:80
ncd_t::l2_z
double l2_z[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Altitude [km].
Definition: diff_apr.c:101
ncd_t::l1_sat_lat
double l1_sat_lat[L1_NTRACK]
Satellite latitude [deg].
Definition: diff_apr.c:92
ncd_t::ncid
int ncid
NetCDF file ID.
Definition: diff_apr.c:71
ncd_t::l1_sat_lon
double l1_sat_lon[L1_NTRACK]
Satellite longitude [deg].
Definition: diff_apr.c:89
ncd_t::l2_p
double l2_p[L2_NLAY]
Pressure [hPa].
Definition: diff_apr.c:104
ncd_t::l2_t
double l2_t[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Temperature [K].
Definition: diff_apr.c:107
ncd_t::l1_nu
double l1_nu[L1_NCHAN]
Channel frequencies [cm^-1].
Definition: diff_apr.c:95
ncd_t::l1_sat_z
double l1_sat_z[L1_NTRACK]
Satellite altitude [km].
Definition: diff_apr.c:86
ncd_t::l1_rad
float l1_rad[L1_NTRACK][L1_NXTRACK][L1_NCHAN]
Radiance [W/(m^2 sr cm^-1)].
Definition: diff_apr.c:98
ncd_t::l1_time
double l1_time[L1_NTRACK][L1_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: diff_apr.c:77
Generated by doxygen 1.9.4