AIRS Code Collection
Data Structures | Macros | Functions
nlte.c File Reference

Retrieval of non-LTE index. More...

#include <mpi.h>
#include <omp.h>
#include <netcdf.h>
#include "jurassic.h"

Go to the source code of this file.

Data Structures

struct  ncd_t
 Buffer for netCDF data. More...
 

Macros

#define NC(cmd)
 Execute netCDF library command and check result. More...
 
#define L1_NCHAN   34
 Number of AIRS radiance channels (don't change). More...
 
#define L1_NTRACK   135
 Along-track size of AIRS radiance granule (don't change). More...
 
#define L1_NXTRACK   90
 Across-track size of AIRS radiance granule (don't change). More...
 
#define L2_NLAY   27
 Number of AIRS pressure layers (don't change). More...
 
#define L2_NTRACK   45
 Along-track size of AIRS retrieval granule (don't change). More...
 
#define L2_NXTRACK   30
 Across-track size of AIRS retrieval granule (don't change). More...
 

Functions

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. More...
 
void fill_gaps (double x[L2_NTRACK][L2_NXTRACK][L2_NLAY], double cx, double cy)
 Fill data gaps in L2 data. More...
 
void init_l2 (ncd_t *ncd, int track, int xtrack, ctl_t *ctl, atm_t *atm)
 Initialize with AIRS Level-2 data. More...
 
void read_nc (char *filename, ncd_t *ncd)
 Read netCDF file. More...
 
int main (int argc, char *argv[])
 

Detailed Description

Retrieval of non-LTE index.

Definition in file nlte.c.

Macro Definition Documentation

◆ NC

#define NC (   cmd)
Value:
{ \
int nc_result=(cmd); \
if(nc_result!=NC_NOERR) \
ERRMSG("%s", nc_strerror(nc_result)); \
}

Execute netCDF library command and check result.

Definition at line 36 of file nlte.c.

◆ L1_NCHAN

#define L1_NCHAN   34

Number of AIRS radiance channels (don't change).

Definition at line 47 of file nlte.c.

◆ L1_NTRACK

#define L1_NTRACK   135

Along-track size of AIRS radiance granule (don't change).

Definition at line 50 of file nlte.c.

◆ L1_NXTRACK

#define L1_NXTRACK   90

Across-track size of AIRS radiance granule (don't change).

Definition at line 53 of file nlte.c.

◆ L2_NLAY

#define L2_NLAY   27

Number of AIRS pressure layers (don't change).

Definition at line 56 of file nlte.c.

◆ L2_NTRACK

#define L2_NTRACK   45

Along-track size of AIRS retrieval granule (don't change).

Definition at line 59 of file nlte.c.

◆ L2_NXTRACK

#define L2_NXTRACK   30

Across-track size of AIRS retrieval granule (don't change).

Definition at line 62 of file nlte.c.

Function Documentation

◆ 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.

Add variable to netCDF file.

Definition at line 403 of file nlte.c.

411 {
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}
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 at line 430 of file nlte.c.

433 {
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}
L2_NXTRACK
#define L2_NXTRACK
Across-track size of AIRS retrieval granule (don't change).
Definition: nlte.c:62
L2_NLAY
#define L2_NLAY
Number of AIRS pressure layers (don't change).
Definition: nlte.c:56
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 at line 475 of file nlte.c.

480 {
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}
ncd_t::l2_z
double l2_z[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Altitude [km].
Definition: diff_apr.c:101
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

◆ read_nc()

void read_nc ( char *  filename,
ncd_t *  ncd 
)

Read netCDF file.

Definition at line 532 of file nlte.c.

534 {
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}
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::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::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

◆ main()

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

Definition at line 159 of file nlte.c.

161 {
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}
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
read_nc
void read_nc(char *filename, ncd_t *ncd)
Read netCDF file.
Definition: nlte.c:532
L1_NCHAN
#define L1_NCHAN
Number of AIRS radiance channels (don't change).
Definition: nlte.c:47
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
Here is the call graph for this function:
Generated by doxygen 1.9.4