IASI Code Collection
extract.c
Go to the documentation of this file.
1/*
2 This file is part of the IASI Code Collection.
3
4 the IASI 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 IASI 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 IASI Code Collection. If not, see
16 <http://www.gnu.org/licenses/>.
17
18 Copyright (C) 2019-2026 Forschungszentrum Juelich GmbH
19*/
20
26#include "libiasi.h"
27
28/* ------------------------------------------------------------
29 Dimensions...
30 ------------------------------------------------------------ */
31
33#define EP 72
34
36#define EX 362
37
39#define EY 182
40
41/* ------------------------------------------------------------
42 Global variables...
43 ------------------------------------------------------------ */
44
46int iasi_chan[L1_NCHAN] = { 71, 88, 89, 90, 91, 92, 93, 94, 95, 96,
47 97, 98, 99, 6712, 6720, 6735, 6742, 6749, 6750,
48 6756, 6757, 6763, 6764, 6770, 6771, 6777, 6778,
49 6784, 6791, 6797, 6838, 6855, 6866
50};
51
52/* ------------------------------------------------------------
53 Structs...
54 ------------------------------------------------------------ */
55
57typedef struct {
58
60 double dt_met;
61
63 char met_geopot[LEN];
64
66 int met_dx;
67
69 int met_dy;
70
72 int met_dp;
73
75 int met_sx;
76
78 int met_sy;
79
81 int met_sp;
82
83} ctl2_t;
84
86typedef struct {
87
89 double time;
90
92 int nx;
93
95 int ny;
96
98 int np;
99
101 double lon[EX];
102
104 double lat[EY];
105
107 double p[EP];
108
110 double ps[EX][EY];
111
113 double pt[EX][EY];
114
116 float z[EX][EY][EP];
117
119 float t[EX][EY][EP];
120
122 float u[EX][EY][EP];
123
125 float v[EX][EY][EP];
126
128 float w[EX][EY][EP];
129
131 float pv[EX][EY][EP];
132
134 float h2o[EX][EY][EP];
135
137 float o3[EX][EY][EP];
138
140 float pl[EX][EY][EP];
141
142} met_t;
143
144/* ------------------------------------------------------------
145 Functions...
146 ------------------------------------------------------------ */
147
149void get_met(
150 ctl2_t * ctl2,
151 char *metbase,
152 double t,
153 met_t * met0,
154 met_t * met1);
155
157void get_met_help(
158 double t,
159 int direct,
160 char *metbase,
161 double dt_met,
162 char *filename);
163
165void intpol_met_2d(
166 double array[EX][EY],
167 int ix,
168 int iy,
169 double wx,
170 double wy,
171 double *var);
172
174void intpol_met_3d(
175 float array[EX][EY][EP],
176 int ip,
177 int ix,
178 int iy,
179 double wp,
180 double wx,
181 double wy,
182 double *var);
183
185void intpol_met_space(
186 met_t * met,
187 double p,
188 double lon,
189 double lat,
190 double *ps,
191 double *pt,
192 double *z,
193 double *t,
194 double *u,
195 double *v,
196 double *w,
197 double *pv,
198 double *h2o,
199 double *o3);
200
202void intpol_met_time(
203 met_t * met0,
204 met_t * met1,
205 double ts,
206 double p,
207 double lon,
208 double lat,
209 double *ps,
210 double *pt,
211 double *z,
212 double *t,
213 double *u,
214 double *v,
215 double *w,
216 double *pv,
217 double *h2o,
218 double *o3);
219
221void read_ctl2(
222 int argc,
223 char *argv[],
224 ctl2_t * ctl2);
225
227void read_met(
228 ctl2_t * ctl2,
229 char *filename,
230 met_t * met);
231
233void read_met_extrapolate(
234 met_t * met);
235
237void read_met_geopot(
238 ctl2_t * ctl2,
239 met_t * met);
240
242void read_met_help(
243 int ncid,
244 char *varname,
245 char *varname2,
246 met_t * met,
247 float dest[EX][EY][EP],
248 float scl);
249
251void read_met_periodic(
252 met_t * met);
253
255void read_met_sample(
256 ctl2_t * ctl2,
257 met_t * met);
258
259/* ------------------------------------------------------------
260 Main...
261 ------------------------------------------------------------ */
262
263int main(
264 int argc,
265 char *argv[]) {
266
267 static iasi_rad_t *iasi_rad;
268
269 static iasi_l1_t l1;
270 static iasi_l2_t l2;
271
272 static ctl2_t ctl2;
273
274 met_t *met0, *met1;
275
276 double ts;
277
278 int ichan, lay, track = 0, xtrack, format;
279
280 /* Check arguments... */
281 if (argc < 4)
282 ERRMSG("Give parameters: <ctl> <iasi_l1_file> <metbase> <out.nc>");
283
284 /* Allocate... */
285 ALLOC(iasi_rad, iasi_rad_t, 1);
286 ALLOC(met0, met_t, 1);
287 ALLOC(met1, met_t, 1);
288
289 /* Read control parameters... */
290 read_ctl2(argc, argv, &ctl2);
291 format = (int) scan_ctl(argc, argv, "FORMAT", -1, "1", NULL);
292
293 /* Read IASI data... */
294 iasi_read(format, argv[2], iasi_rad);
295
296 /* Copy data to struct... */
297 l1.ntrack = (size_t) iasi_rad->ntrack;
298 for (track = 0; track < iasi_rad->ntrack; track++)
299 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
300 l1.time[track][xtrack]
301 = iasi_rad->Time[track][xtrack];
302 l1.lon[track][xtrack]
303 = iasi_rad->Longitude[track][xtrack];
304 l1.lat[track][xtrack]
305 = iasi_rad->Latitude[track][xtrack];
306 l1.sat_z[track]
307 = iasi_rad->Sat_z[track];
308 l1.sat_lon[track]
309 = iasi_rad->Sat_lon[track];
310 l1.sat_lat[track]
311 = iasi_rad->Sat_lat[track];
312 for (ichan = 0; ichan < L1_NCHAN; ichan++) {
313 l1.nu[ichan]
314 = iasi_rad->freq[iasi_chan[ichan]];
315 l1.rad[track][xtrack][ichan]
316 = iasi_rad->Rad[track][xtrack][iasi_chan[ichan]];
317 }
318 }
319
320 /* Write netCDF file... */
321 write_l1(argv[4], &l1);
322
323 /* Read meteo data... */
324 l2.ntrack = l1.ntrack;
325 ts = l1.time[l2.ntrack / 2][L1_NXTRACK / 2];
326 get_met(&ctl2, argv[3], ts, met0, met1);
327
328 /* Interpolate meteo data... */
329 for (track = 0; track < (int) l2.ntrack; track++)
330 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++)
331 for (lay = 0; lay < L2_NLAY; lay++) {
332 l2.time[track][xtrack] = l1.time[track][xtrack];
333 l2.lon[track][xtrack] = l1.lon[track][xtrack];
334 l2.lat[track][xtrack] = l1.lat[track][xtrack];
335 l2.p[lay] = 1013.25 * exp(-2.5 * lay / 7.0);
336 intpol_met_time(met0, met1, ts, l2.p[lay],
337 l2.lon[track][xtrack], l2.lat[track][xtrack],
338 NULL, NULL, &l2.z[track][xtrack][lay],
339 &l2.t[track][xtrack][lay],
340 NULL, NULL, NULL, NULL, NULL, NULL);
341 }
342
343 /* Write netCDF file... */
344 write_l2(argv[4], &l2);
345
346 /* Free... */
347 free(iasi_rad);
348 free(met0);
349 free(met1);
350
351 return EXIT_SUCCESS;
352}
353
354/*****************************************************************************/
355
356void get_met(
357 ctl2_t *ctl2,
358 char *metbase,
359 double t,
360 met_t *met0,
361 met_t *met1) {
362
363 char filename[LEN];
364
365 static int init;
366
367 /* Init... */
368 if (!init) {
369 init = 1;
370
371 get_met_help(t, -1, metbase, ctl2->dt_met, filename);
372 read_met(ctl2, filename, met0);
373
374 get_met_help(t + 1.0, 1, metbase, ctl2->dt_met, filename);
375 read_met(ctl2, filename, met1);
376 }
377
378 /* Read new data for forward trajectories... */
379 if (t > met1->time) {
380 memcpy(met0, met1, sizeof(met_t));
381 get_met_help(t, 1, metbase, ctl2->dt_met, filename);
382 read_met(ctl2, filename, met1);
383 }
384}
385
386/*****************************************************************************/
387
388void get_met_help(
389 double t,
390 int direct,
391 char *metbase,
392 double dt_met,
393 char *filename) {
394
395 double t6, r;
396
397 int year, mon, day, hour, min, sec;
398
399 /* Round time to fixed intervals... */
400 if (direct == -1)
401 t6 = floor(t / dt_met) * dt_met;
402 else
403 t6 = ceil(t / dt_met) * dt_met;
404
405 /* Decode time... */
406 jsec2time(t6, &year, &mon, &day, &hour, &min, &sec, &r);
407
408 /* Set filename... */
409 sprintf(filename, "%s_%d_%02d_%02d_%02d.nc", metbase, year, mon, day, hour);
410}
411
412/*****************************************************************************/
413
414void intpol_met_2d(
415 double array[EX][EY],
416 int ix,
417 int iy,
418 double wx,
419 double wy,
420 double *var) {
421
422 double aux00, aux01, aux10, aux11;
423
424 /* Set variables... */
425 aux00 = array[ix][iy];
426 aux01 = array[ix][iy + 1];
427 aux10 = array[ix + 1][iy];
428 aux11 = array[ix + 1][iy + 1];
429
430 /* Interpolate horizontally... */
431 aux00 = wy * (aux00 - aux01) + aux01;
432 aux11 = wy * (aux10 - aux11) + aux11;
433 *var = wx * (aux00 - aux11) + aux11;
434}
435
436/*****************************************************************************/
437
438void intpol_met_3d(
439 float array[EX][EY][EP],
440 int ip,
441 int ix,
442 int iy,
443 double wp,
444 double wx,
445 double wy,
446 double *var) {
447
448 double aux00, aux01, aux10, aux11;
449
450 /* Interpolate vertically... */
451 aux00 = wp * (array[ix][iy][ip] - array[ix][iy][ip + 1])
452 + array[ix][iy][ip + 1];
453 aux01 = wp * (array[ix][iy + 1][ip] - array[ix][iy + 1][ip + 1])
454 + array[ix][iy + 1][ip + 1];
455 aux10 = wp * (array[ix + 1][iy][ip] - array[ix + 1][iy][ip + 1])
456 + array[ix + 1][iy][ip + 1];
457 aux11 = wp * (array[ix + 1][iy + 1][ip] - array[ix + 1][iy + 1][ip + 1])
458 + array[ix + 1][iy + 1][ip + 1];
459
460 /* Interpolate horizontally... */
461 aux00 = wy * (aux00 - aux01) + aux01;
462 aux11 = wy * (aux10 - aux11) + aux11;
463 *var = wx * (aux00 - aux11) + aux11;
464}
465
466/*****************************************************************************/
467
468void intpol_met_space(
469 met_t *met,
470 double p,
471 double lon,
472 double lat,
473 double *ps,
474 double *pt,
475 double *z,
476 double *t,
477 double *u,
478 double *v,
479 double *w,
480 double *pv,
481 double *h2o,
482 double *o3) {
483
484 double wp, wx, wy;
485
486 int ip, ix, iy;
487
488 /* Check longitude... */
489 if (met->lon[met->nx - 1] > 180 && lon < 0)
490 lon += 360;
491
492 /* Get indices... */
493 ip = locate_irr(met->p, met->np, p);
494 ix = locate_reg(met->lon, met->nx, lon);
495 iy = locate_reg(met->lat, met->ny, lat);
496
497 /* Get weights... */
498 wp = (met->p[ip + 1] - p) / (met->p[ip + 1] - met->p[ip]);
499 wx = (met->lon[ix + 1] - lon) / (met->lon[ix + 1] - met->lon[ix]);
500 wy = (met->lat[iy + 1] - lat) / (met->lat[iy + 1] - met->lat[iy]);
501
502 /* Interpolate... */
503 if (ps != NULL)
504 intpol_met_2d(met->ps, ix, iy, wx, wy, ps);
505 if (pt != NULL)
506 intpol_met_2d(met->pt, ix, iy, wx, wy, pt);
507 if (z != NULL)
508 intpol_met_3d(met->z, ip, ix, iy, wp, wx, wy, z);
509 if (t != NULL)
510 intpol_met_3d(met->t, ip, ix, iy, wp, wx, wy, t);
511 if (u != NULL)
512 intpol_met_3d(met->u, ip, ix, iy, wp, wx, wy, u);
513 if (v != NULL)
514 intpol_met_3d(met->v, ip, ix, iy, wp, wx, wy, v);
515 if (w != NULL)
516 intpol_met_3d(met->w, ip, ix, iy, wp, wx, wy, w);
517 if (pv != NULL)
518 intpol_met_3d(met->pv, ip, ix, iy, wp, wx, wy, pv);
519 if (h2o != NULL)
520 intpol_met_3d(met->h2o, ip, ix, iy, wp, wx, wy, h2o);
521 if (o3 != NULL)
522 intpol_met_3d(met->o3, ip, ix, iy, wp, wx, wy, o3);
523}
524
525/*****************************************************************************/
526
527void intpol_met_time(
528 met_t *met0,
529 met_t *met1,
530 double ts,
531 double p,
532 double lon,
533 double lat,
534 double *ps,
535 double *pt,
536 double *z,
537 double *t,
538 double *u,
539 double *v,
540 double *w,
541 double *pv,
542 double *h2o,
543 double *o3) {
544
545 double h2o0, h2o1, o30, o31, ps0, ps1, pt0, pt1, pv0, pv1, t0, t1, u0, u1,
546 v0, v1, w0, w1, wt, z0, z1;
547
548 /* Spatial interpolation... */
549 intpol_met_space(met0, p, lon, lat,
550 ps == NULL ? NULL : &ps0,
551 pt == NULL ? NULL : &pt0,
552 z == NULL ? NULL : &z0,
553 t == NULL ? NULL : &t0,
554 u == NULL ? NULL : &u0,
555 v == NULL ? NULL : &v0,
556 w == NULL ? NULL : &w0,
557 pv == NULL ? NULL : &pv0,
558 h2o == NULL ? NULL : &h2o0, o3 == NULL ? NULL : &o30);
559 intpol_met_space(met1, p, lon, lat,
560 ps == NULL ? NULL : &ps1,
561 pt == NULL ? NULL : &pt1,
562 z == NULL ? NULL : &z1,
563 t == NULL ? NULL : &t1,
564 u == NULL ? NULL : &u1,
565 v == NULL ? NULL : &v1,
566 w == NULL ? NULL : &w1,
567 pv == NULL ? NULL : &pv1,
568 h2o == NULL ? NULL : &h2o1, o3 == NULL ? NULL : &o31);
569
570 /* Get weighting factor... */
571 wt = (met1->time - ts) / (met1->time - met0->time);
572
573 /* Interpolate... */
574 if (ps != NULL)
575 *ps = wt * (ps0 - ps1) + ps1;
576 if (pt != NULL)
577 *pt = wt * (pt0 - pt1) + pt1;
578 if (z != NULL)
579 *z = wt * (z0 - z1) + z1;
580 if (t != NULL)
581 *t = wt * (t0 - t1) + t1;
582 if (u != NULL)
583 *u = wt * (u0 - u1) + u1;
584 if (v != NULL)
585 *v = wt * (v0 - v1) + v1;
586 if (w != NULL)
587 *w = wt * (w0 - w1) + w1;
588 if (pv != NULL)
589 *pv = wt * (pv0 - pv1) + pv1;
590 if (h2o != NULL)
591 *h2o = wt * (h2o0 - h2o1) + h2o1;
592 if (o3 != NULL)
593 *o3 = wt * (o30 - o31) + o31;
594}
595
596/*****************************************************************************/
597
598void read_ctl2(
599 int argc,
600 char *argv[],
601 ctl2_t *ctl2) {
602
603 /* Meteorological data... */
604 ctl2->dt_met = scan_ctl(argc, argv, "DT_MET", -1, "21600", NULL);
605 scan_ctl(argc, argv, "MET_GEOPOT", -1, "", ctl2->met_geopot);
606 ctl2->met_dx = (int) scan_ctl(argc, argv, "MET_DX", -1, "1", NULL);
607 ctl2->met_dy = (int) scan_ctl(argc, argv, "MET_DY", -1, "1", NULL);
608 ctl2->met_dp = (int) scan_ctl(argc, argv, "MET_DP", -1, "1", NULL);
609 ctl2->met_sx = (int) scan_ctl(argc, argv, "MET_SX", -1, "20", NULL);
610 ctl2->met_sy = (int) scan_ctl(argc, argv, "MET_SY", -1, "10", NULL);
611 ctl2->met_sp = (int) scan_ctl(argc, argv, "MET_SP", -1, "1", NULL);
612}
613
614/*****************************************************************************/
615
616void read_met(
617 ctl2_t *ctl2,
618 char *filename,
619 met_t *met) {
620
621 char levname[LEN], tstr[10];
622
623 static float help[EX * EY];
624
625 int ix, iy, ip, dimid, ncid, varid, year, mon, day, hour;
626
627 size_t np, nx, ny;
628
629 /* Write info... */
630 printf("Read meteorological data: %s\n", filename);
631
632 /* Get time from filename... */
633 sprintf(tstr, "%.4s", &filename[strlen(filename) - 16]);
634 year = atoi(tstr);
635 sprintf(tstr, "%.2s", &filename[strlen(filename) - 11]);
636 mon = atoi(tstr);
637 sprintf(tstr, "%.2s", &filename[strlen(filename) - 8]);
638 day = atoi(tstr);
639 sprintf(tstr, "%.2s", &filename[strlen(filename) - 5]);
640 hour = atoi(tstr);
641 time2jsec(year, mon, day, hour, 0, 0, 0, &met->time);
642
643 /* Open netCDF file... */
644 NC(nc_open(filename, NC_NOWRITE, &ncid));
645
646 /* Get dimensions... */
647 NC(nc_inq_dimid(ncid, "lon", &dimid));
648 NC(nc_inq_dimlen(ncid, dimid, &nx));
649 if (nx < 2 || nx > EX)
650 ERRMSG("Number of longitudes out of range!");
651
652 NC(nc_inq_dimid(ncid, "lat", &dimid));
653 NC(nc_inq_dimlen(ncid, dimid, &ny));
654 if (ny < 2 || ny > EY)
655 ERRMSG("Number of latitudes out of range!");
656
657 sprintf(levname, "lev");
658 NC(nc_inq_dimid(ncid, levname, &dimid));
659 NC(nc_inq_dimlen(ncid, dimid, &np));
660 if (np == 1) {
661 sprintf(levname, "lev_2");
662 NC(nc_inq_dimid(ncid, levname, &dimid));
663 NC(nc_inq_dimlen(ncid, dimid, &np));
664 }
665 if (np < 2 || np > EP)
666 ERRMSG("Number of levels out of range!");
667
668 /* Store dimensions... */
669 met->np = (int) np;
670 met->nx = (int) nx;
671 met->ny = (int) ny;
672
673 /* Get horizontal grid... */
674 NC(nc_inq_varid(ncid, "lon", &varid));
675 NC(nc_get_var_double(ncid, varid, met->lon));
676 NC(nc_inq_varid(ncid, "lat", &varid));
677 NC(nc_get_var_double(ncid, varid, met->lat));
678
679 /* Read meteorological data... */
680 read_met_help(ncid, "t", "T", met, met->t, 1.0);
681 read_met_help(ncid, "u", "U", met, met->u, 1.0);
682 read_met_help(ncid, "v", "V", met, met->v, 1.0);
683 read_met_help(ncid, "w", "W", met, met->w, 0.01f);
684 read_met_help(ncid, "q", "Q", met, met->h2o, 1.608f);
685 read_met_help(ncid, "o3", "O3", met, met->o3, 0.602f);
686
687 /* Read pressure levels from file... */
688 NC(nc_inq_varid(ncid, levname, &varid));
689 NC(nc_get_var_double(ncid, varid, met->p));
690 for (ip = 0; ip < met->np; ip++)
691 met->p[ip] /= 100.;
692
693 /* Extrapolate data for lower boundary... */
694 read_met_extrapolate(met);
695
696 /* Check ordering of pressure levels... */
697 for (ip = 1; ip < met->np; ip++)
698 if (met->p[ip - 1] < met->p[ip])
699 ERRMSG("Pressure levels must be descending!");
700
701 /* Read surface pressure... */
702 if (nc_inq_varid(ncid, "ps", &varid) == NC_NOERR
703 || nc_inq_varid(ncid, "PS", &varid) == NC_NOERR) {
704 NC(nc_get_var_float(ncid, varid, help));
705 for (iy = 0; iy < met->ny; iy++)
706 for (ix = 0; ix < met->nx; ix++)
707 met->ps[ix][iy] = help[iy * met->nx + ix] / 100.;
708 } else if (nc_inq_varid(ncid, "lnsp", &varid) == NC_NOERR
709 || nc_inq_varid(ncid, "LNSP", &varid) == NC_NOERR) {
710 NC(nc_get_var_float(ncid, varid, help));
711 for (iy = 0; iy < met->ny; iy++)
712 for (ix = 0; ix < met->nx; ix++)
713 met->ps[ix][iy] = exp(help[iy * met->nx + ix]) / 100.;
714 } else
715 for (ix = 0; ix < met->nx; ix++)
716 for (iy = 0; iy < met->ny; iy++)
717 met->ps[ix][iy] = met->p[0];
718
719 /* Create periodic boundary conditions... */
720 read_met_periodic(met);
721
722 /* Calculate geopotential heights... */
723 read_met_geopot(ctl2, met);
724
725 /* Downsampling... */
726 read_met_sample(ctl2, met);
727
728 /* Close file... */
729 NC(nc_close(ncid));
730}
731
732/*****************************************************************************/
733
734void read_met_extrapolate(
735 met_t *met) {
736
737 int ip, ip0, ix, iy;
738
739 /* Loop over columns... */
740 for (ix = 0; ix < met->nx; ix++)
741 for (iy = 0; iy < met->ny; iy++) {
742
743 /* Find lowest valid data point... */
744 for (ip0 = met->np - 1; ip0 >= 0; ip0--)
745 if (!gsl_finite(met->t[ix][iy][ip0])
746 || !gsl_finite(met->u[ix][iy][ip0])
747 || !gsl_finite(met->v[ix][iy][ip0])
748 || !gsl_finite(met->w[ix][iy][ip0]))
749 break;
750
751 /* Extrapolate... */
752 for (ip = ip0; ip >= 0; ip--) {
753 met->t[ix][iy][ip] = met->t[ix][iy][ip + 1];
754 met->u[ix][iy][ip] = met->u[ix][iy][ip + 1];
755 met->v[ix][iy][ip] = met->v[ix][iy][ip + 1];
756 met->w[ix][iy][ip] = met->w[ix][iy][ip + 1];
757 met->h2o[ix][iy][ip] = met->h2o[ix][iy][ip + 1];
758 met->o3[ix][iy][ip] = met->o3[ix][iy][ip + 1];
759 }
760 }
761}
762
763/*****************************************************************************/
764
765void read_met_geopot(
766 ctl2_t *ctl2,
767 met_t *met) {
768
769 static double topo_lat[EY], topo_lon[EX], topo_z[EX][EY];
770
771 static int init, topo_nx = -1, topo_ny;
772
773 FILE *in;
774
775 char line[LEN];
776
777 double data[30], lat, lon, rlat, rlon, rlon_old = -999, rz, ts, z0, z1;
778
779 float help[EX][EY];
780
781 int ip, ip0, ix, ix2, ix3, iy, iy2, n, tx, ty;
782
783 /* Initialize geopotential heights... */
784 for (ix = 0; ix < met->nx; ix++)
785 for (iy = 0; iy < met->ny; iy++)
786 for (ip = 0; ip < met->np; ip++)
787 met->z[ix][iy][ip] = GSL_NAN;
788
789 /* Check filename... */
790 if (ctl2->met_geopot[0] == '-')
791 return;
792
793 /* Read surface geopotential... */
794 if (!init) {
795
796 /* Write info... */
797 printf("Read surface geopotential: %s\n", ctl2->met_geopot);
798
799 /* Open file... */
800 if (!(in = fopen(ctl2->met_geopot, "r")))
801 ERRMSG("Cannot open file!");
802
803 /* Read data... */
804 while (fgets(line, LEN, in))
805 if (sscanf(line, "%lg %lg %lg", &rlon, &rlat, &rz) == 3) {
806 if (rlon != rlon_old) {
807 if ((++topo_nx) >= EX)
808 ERRMSG("Too many longitudes!");
809 topo_ny = 0;
810 }
811 rlon_old = rlon;
812 topo_lon[topo_nx] = rlon;
813 topo_lat[topo_ny] = rlat;
814 topo_z[topo_nx][topo_ny] = rz;
815 if ((++topo_ny) >= EY)
816 ERRMSG("Too many latitudes!");
817 }
818 if ((++topo_nx) >= EX)
819 ERRMSG("Too many longitudes!");
820
821 /* Close file... */
822 fclose(in);
823
824 /* Check grid spacing... */
825 if (fabs(met->lon[0] - met->lon[1]) != fabs(topo_lon[0] - topo_lon[1])
826 || fabs(met->lat[0] - met->lat[1]) != fabs(topo_lat[0] - topo_lat[1]))
827 printf("Warning: Grid spacing does not match!\n");
828
829 /* Set init flag... */
830 init = 1;
831 }
832
833 /* Apply hydrostatic equation to calculate geopotential heights... */
834 for (ix = 0; ix < met->nx; ix++)
835 for (iy = 0; iy < met->ny; iy++) {
836
837 /* Get surface height... */
838 lon = met->lon[ix];
839 if (lon < topo_lon[0])
840 lon += 360;
841 else if (lon > topo_lon[topo_nx - 1])
842 lon -= 360;
843 lat = met->lat[iy];
844 tx = locate_reg(topo_lon, topo_nx, lon);
845 ty = locate_reg(topo_lat, topo_ny, lat);
846 z0 = LIN(topo_lon[tx], topo_z[tx][ty],
847 topo_lon[tx + 1], topo_z[tx + 1][ty], lon);
848 z1 = LIN(topo_lon[tx], topo_z[tx][ty + 1],
849 topo_lon[tx + 1], topo_z[tx + 1][ty + 1], lon);
850 z0 = LIN(topo_lat[ty], z0, topo_lat[ty + 1], z1, lat);
851
852 /* Find surface pressure level... */
853 ip0 = locate_irr(met->p, met->np, met->ps[ix][iy]);
854
855 /* Get surface temperature... */
856 ts = LIN(met->p[ip0], met->t[ix][iy][ip0],
857 met->p[ip0 + 1], met->t[ix][iy][ip0 + 1], met->ps[ix][iy]);
858
859 /* Upper part of profile... */
860 met->z[ix][iy][ip0 + 1]
861 = (float) (z0 + 8.31441 / 28.9647 / G0
862 * 0.5 * (ts + met->t[ix][iy][ip0 + 1])
863 * log(met->ps[ix][iy] / met->p[ip0 + 1]));
864 for (ip = ip0 + 2; ip < met->np; ip++)
865 met->z[ix][iy][ip]
866 = (float) (met->z[ix][iy][ip - 1] + 8.31441 / 28.9647 / G0
867 * 0.5 * (met->t[ix][iy][ip - 1] + met->t[ix][iy][ip])
868 * log(met->p[ip - 1] / met->p[ip]));
869 }
870
871 /* Smooth fields... */
872 for (ip = 0; ip < met->np; ip++) {
873
874 /* Median filter... */
875 for (ix = 0; ix < met->nx; ix++)
876 for (iy = 0; iy < met->nx; iy++) {
877 n = 0;
878 for (ix2 = ix - 2; ix2 <= ix + 2; ix2++) {
879 ix3 = ix2;
880 if (ix3 < 0)
881 ix3 += met->nx;
882 if (ix3 >= met->nx)
883 ix3 -= met->nx;
884 for (iy2 = GSL_MAX(iy - 2, 0); iy2 <= GSL_MIN(iy + 2, met->ny - 1);
885 iy2++)
886 if (gsl_finite(met->z[ix3][iy2][ip])) {
887 data[n] = met->z[ix3][iy2][ip];
888 n++;
889 }
890 }
891 if (n > 0) {
892 gsl_sort(data, 1, (size_t) n);
893 help[ix][iy] = (float)
894 gsl_stats_median_from_sorted_data(data, 1, (size_t) n);
895 } else
896 help[ix][iy] = GSL_NAN;
897 }
898
899 /* Copy data... */
900 for (ix = 0; ix < met->nx; ix++)
901 for (iy = 0; iy < met->nx; iy++)
902 met->z[ix][iy][ip] = help[ix][iy];
903 }
904}
905
906/*****************************************************************************/
907
908void read_met_help(
909 int ncid,
910 char *varname,
911 char *varname2,
912 met_t *met,
913 float dest[EX][EY][EP],
914 float scl) {
915
916 static float help[EX * EY * EP];
917
918 int ip, ix, iy, varid;
919
920 /* Check if variable exists... */
921 if (nc_inq_varid(ncid, varname, &varid) != NC_NOERR)
922 if (nc_inq_varid(ncid, varname2, &varid) != NC_NOERR)
923 return;
924
925 /* Read data... */
926 NC(nc_get_var_float(ncid, varid, help));
927
928 /* Copy and check data... */
929 for (ix = 0; ix < met->nx; ix++)
930 for (iy = 0; iy < met->ny; iy++)
931 for (ip = 0; ip < met->np; ip++) {
932 dest[ix][iy][ip] = help[(ip * met->ny + iy) * met->nx + ix];
933 if (fabsf(dest[ix][iy][ip]) < 1e14f)
934 dest[ix][iy][ip] *= scl;
935 else
936 dest[ix][iy][ip] = GSL_NAN;
937 }
938}
939
940/*****************************************************************************/
941
942void read_met_periodic(
943 met_t *met) {
944
945 int ip, iy;
946
947 /* Check longitudes... */
948 if (!(fabs(met->lon[met->nx - 1] - met->lon[0]
949 + met->lon[1] - met->lon[0] - 360) < 0.01))
950 return;
951
952 /* Increase longitude counter... */
953 if ((++met->nx) > EX)
954 ERRMSG("Cannot create periodic boundary conditions!");
955
956 /* Set longitude... */
957 met->lon[met->nx - 1] = met->lon[met->nx - 2] + met->lon[1] - met->lon[0];
958
959 /* Loop over latitudes and pressure levels... */
960 for (iy = 0; iy < met->ny; iy++)
961 for (ip = 0; ip < met->np; ip++) {
962 met->ps[met->nx - 1][iy] = met->ps[0][iy];
963 met->pt[met->nx - 1][iy] = met->pt[0][iy];
964 met->z[met->nx - 1][iy][ip] = met->z[0][iy][ip];
965 met->t[met->nx - 1][iy][ip] = met->t[0][iy][ip];
966 met->u[met->nx - 1][iy][ip] = met->u[0][iy][ip];
967 met->v[met->nx - 1][iy][ip] = met->v[0][iy][ip];
968 met->w[met->nx - 1][iy][ip] = met->w[0][iy][ip];
969 met->pv[met->nx - 1][iy][ip] = met->pv[0][iy][ip];
970 met->h2o[met->nx - 1][iy][ip] = met->h2o[0][iy][ip];
971 met->o3[met->nx - 1][iy][ip] = met->o3[0][iy][ip];
972 }
973}
974
975/*****************************************************************************/
976
977void read_met_sample(
978 ctl2_t *ctl2,
979 met_t *met) {
980
981 met_t *help;
982
983 float w, wsum;
984
985 int ip, ip2, ix, ix2, ix3, iy, iy2;
986
987 /* Check parameters... */
988 if (ctl2->met_dp <= 1 && ctl2->met_dx <= 1 && ctl2->met_dy <= 1
989 && ctl2->met_sp <= 1 && ctl2->met_sx <= 1 && ctl2->met_sy <= 1)
990 return;
991
992 /* Allocate... */
993 ALLOC(help, met_t, 1);
994
995 /* Copy data... */
996 help->nx = met->nx;
997 help->ny = met->ny;
998 help->np = met->np;
999 memcpy(help->lon, met->lon, sizeof(met->lon));
1000 memcpy(help->lat, met->lat, sizeof(met->lat));
1001 memcpy(help->p, met->p, sizeof(met->p));
1002
1003 /* Smoothing... */
1004 for (ix = 0; ix < met->nx; ix += ctl2->met_dx) {
1005 for (iy = 0; iy < met->ny; iy += ctl2->met_dy) {
1006 for (ip = 0; ip < met->np; ip += ctl2->met_dp) {
1007 help->ps[ix][iy] = 0;
1008 help->pt[ix][iy] = 0;
1009 help->z[ix][iy][ip] = 0;
1010 help->t[ix][iy][ip] = 0;
1011 help->u[ix][iy][ip] = 0;
1012 help->v[ix][iy][ip] = 0;
1013 help->w[ix][iy][ip] = 0;
1014 help->pv[ix][iy][ip] = 0;
1015 help->h2o[ix][iy][ip] = 0;
1016 help->o3[ix][iy][ip] = 0;
1017 wsum = 0;
1018 for (ix2 = ix - ctl2->met_sx + 1; ix2 <= ix + ctl2->met_sx - 1; ix2++) {
1019 ix3 = ix2;
1020 if (ix3 < 0)
1021 ix3 += met->nx;
1022 else if (ix3 >= met->nx)
1023 ix3 -= met->nx;
1024
1025 for (iy2 = GSL_MAX(iy - ctl2->met_sy + 1, 0);
1026 iy2 <= GSL_MIN(iy + ctl2->met_sy - 1, met->ny - 1); iy2++)
1027 for (ip2 = GSL_MAX(ip - ctl2->met_sp + 1, 0);
1028 ip2 <= GSL_MIN(ip + ctl2->met_sp - 1, met->np - 1); ip2++) {
1029 w = (1.0f - (float) abs(ix - ix2) / (float) ctl2->met_sx)
1030 * (1.0f - (float) abs(iy - iy2) / (float) ctl2->met_sy)
1031 * (1.0f - (float) abs(ip - ip2) / (float) ctl2->met_sp);
1032 help->ps[ix][iy] += w * met->ps[ix3][iy2];
1033 help->pt[ix][iy] += w * met->pt[ix3][iy2];
1034 help->z[ix][iy][ip] += w * met->z[ix3][iy2][ip2];
1035 help->t[ix][iy][ip] += w * met->t[ix3][iy2][ip2];
1036 help->u[ix][iy][ip] += w * met->u[ix3][iy2][ip2];
1037 help->v[ix][iy][ip] += w * met->v[ix3][iy2][ip2];
1038 help->w[ix][iy][ip] += w * met->w[ix3][iy2][ip2];
1039 help->pv[ix][iy][ip] += w * met->pv[ix3][iy2][ip2];
1040 help->h2o[ix][iy][ip] += w * met->h2o[ix3][iy2][ip2];
1041 help->o3[ix][iy][ip] += w * met->o3[ix3][iy2][ip2];
1042 wsum += w;
1043 }
1044 }
1045 help->ps[ix][iy] /= wsum;
1046 help->pt[ix][iy] /= wsum;
1047 help->t[ix][iy][ip] /= wsum;
1048 help->z[ix][iy][ip] /= wsum;
1049 help->u[ix][iy][ip] /= wsum;
1050 help->v[ix][iy][ip] /= wsum;
1051 help->w[ix][iy][ip] /= wsum;
1052 help->pv[ix][iy][ip] /= wsum;
1053 help->h2o[ix][iy][ip] /= wsum;
1054 help->o3[ix][iy][ip] /= wsum;
1055 }
1056 }
1057 }
1058
1059 /* Downsampling... */
1060 met->nx = 0;
1061 for (ix = 0; ix < help->nx; ix += ctl2->met_dx) {
1062 met->lon[met->nx] = help->lon[ix];
1063 met->ny = 0;
1064 for (iy = 0; iy < help->ny; iy += ctl2->met_dy) {
1065 met->lat[met->ny] = help->lat[iy];
1066 met->ps[met->nx][met->ny] = help->ps[ix][iy];
1067 met->pt[met->nx][met->ny] = help->pt[ix][iy];
1068 met->np = 0;
1069 for (ip = 0; ip < help->np; ip += ctl2->met_dp) {
1070 met->p[met->np] = help->p[ip];
1071 met->z[met->nx][met->ny][met->np] = help->z[ix][iy][ip];
1072 met->t[met->nx][met->ny][met->np] = help->t[ix][iy][ip];
1073 met->u[met->nx][met->ny][met->np] = help->u[ix][iy][ip];
1074 met->v[met->nx][met->ny][met->np] = help->v[ix][iy][ip];
1075 met->w[met->nx][met->ny][met->np] = help->w[ix][iy][ip];
1076 met->pv[met->nx][met->ny][met->np] = help->pv[ix][iy][ip];
1077 met->h2o[met->nx][met->ny][met->np] = help->h2o[ix][iy][ip];
1078 met->o3[met->nx][met->ny][met->np] = help->o3[ix][iy][ip];
1079 met->np++;
1080 }
1081 met->ny++;
1082 }
1083 met->nx++;
1084 }
1085
1086 /* Free... */
1087 free(help);
1088}
read_ctl2
void read_ctl2(int argc, char *argv[], ctl2_t *ctl2)
Read control parameters.
Definition: extract.c:598
read_met_extrapolate
void read_met_extrapolate(met_t *met)
Extrapolate meteorological data at lower boundary.
Definition: extract.c:734
main
int main(int argc, char *argv[])
Definition: extract.c:263
read_met_periodic
void read_met_periodic(met_t *met)
Create meteorological data with periodic boundary conditions.
Definition: extract.c:942
intpol_met_space
void intpol_met_space(met_t *met, double p, double lon, double lat, double *ps, double *pt, double *z, double *t, double *u, double *v, double *w, double *pv, double *h2o, double *o3)
Spatial interpolation of meteorological data.
Definition: extract.c:468
EY
#define EY
Maximum number of latitudes for meteorological data.
Definition: extract.c:39
intpol_met_2d
void intpol_met_2d(double array[EX][EY], int ix, int iy, double wx, double wy, double *var)
Linear interpolation of 2-D meteorological data.
Definition: extract.c:414
iasi_chan
int iasi_chan[L1_NCHAN]
List of IASI channels (don't change).
Definition: extract.c:46
EX
#define EX
Maximum number of longitudes for meteorological data.
Definition: extract.c:36
intpol_met_3d
void intpol_met_3d(float array[EX][EY][EP], int ip, int ix, int iy, double wp, double wx, double wy, double *var)
Linear interpolation of 3-D meteorological data.
Definition: extract.c:438
read_met_sample
void read_met_sample(ctl2_t *ctl2, met_t *met)
Downsampling of meteorological data.
Definition: extract.c:977
get_met_help
void get_met_help(double t, int direct, char *metbase, double dt_met, char *filename)
Get meteorological data for timestep.
Definition: extract.c:388
read_met_help
void read_met_help(int ncid, char *varname, char *varname2, met_t *met, float dest[EX][EY][EP], float scl)
Read and convert variable from meteorological data file.
Definition: extract.c:908
get_met
void get_met(ctl2_t *ctl2, char *metbase, double t, met_t *met0, met_t *met1)
Get meteorological data for given timestep.
Definition: extract.c:356
EP
#define EP
Maximum number of pressure levels for meteorological data.
Definition: extract.c:33
read_met_geopot
void read_met_geopot(ctl2_t *ctl2, met_t *met)
Calculate geopotential heights.
Definition: extract.c:765
read_met
void read_met(ctl2_t *ctl2, char *filename, met_t *met)
Read meteorological data file.
Definition: extract.c:616
intpol_met_time
void intpol_met_time(met_t *met0, met_t *met1, double ts, double p, double lon, double lat, double *ps, double *pt, double *z, double *t, double *u, double *v, double *w, double *pv, double *h2o, double *o3)
Temporal interpolation of meteorological data.
Definition: extract.c:527
iasi_read
void iasi_read(int format, char *filename, iasi_rad_t *iasi_rad)
Read IASI Level-1 data.
Definition: libiasi.c:297
write_l2
void write_l2(char *filename, iasi_l2_t *l2)
Write IASI Level-2 data.
Definition: libiasi.c:1139
write_l1
void write_l1(char *filename, iasi_l1_t *l1)
Write IASI Level-1 data.
Definition: libiasi.c:1081
libiasi.h
IASI Code Collection library declarations.
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libiasi.h:155
L1_NXTRACK
#define L1_NXTRACK
Across-track size of IASI radiance granule (don't change).
Definition: libiasi.h:102
L2_NLAY
#define L2_NLAY
Number of IASI pressure layers (don't change).
Definition: libiasi.h:105
L1_NCHAN
#define L1_NCHAN
Number of IASI radiance channels (don't change).
Definition: libiasi.h:96
ctl2_t
Control parameters.
Definition: extract.c:57
ctl2_t::met_sy
int met_sy
Smoothing for latitudes.
Definition: extract.c:78
ctl2_t::met_geopot
char met_geopot[LEN]
Surface geopotential data file.
Definition: extract.c:63
ctl2_t::met_sx
int met_sx
Smoothing for longitudes.
Definition: extract.c:75
ctl2_t::met_dp
int met_dp
Stride for pressure levels.
Definition: extract.c:72
ctl2_t::met_sp
int met_sp
Smoothing for pressure levels.
Definition: extract.c:81
ctl2_t::dt_met
double dt_met
Time step of meteorological data [s].
Definition: extract.c:60
ctl2_t::met_dy
int met_dy
Stride for latitudes.
Definition: extract.c:69
ctl2_t::met_dx
int met_dx
Stride for longitudes.
Definition: extract.c:66
iasi_l1_t
IASI Level-1 data.
Definition: libiasi.h:166
iasi_l1_t::sat_z
double sat_z[L1_NTRACK]
Satellite altitude [km].
Definition: libiasi.h:181
iasi_l1_t::nu
double nu[L1_NCHAN]
Channel frequencies [cm^-1].
Definition: libiasi.h:190
iasi_l1_t::sat_lon
double sat_lon[L1_NTRACK]
Satellite longitude [deg].
Definition: libiasi.h:184
iasi_l1_t::sat_lat
double sat_lat[L1_NTRACK]
Satellite latitude [deg].
Definition: libiasi.h:187
iasi_l1_t::lon
double lon[L1_NTRACK][L1_NXTRACK]
Footprint longitude [deg].
Definition: libiasi.h:175
iasi_l1_t::time
double time[L1_NTRACK][L1_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:172
iasi_l1_t::lat
double lat[L1_NTRACK][L1_NXTRACK]
Footprint latitude [deg].
Definition: libiasi.h:178
iasi_l1_t::rad
float rad[L1_NTRACK][L1_NXTRACK][L1_NCHAN]
Radiance [W/(m^2 sr cm^-1)].
Definition: libiasi.h:193
iasi_l1_t::ntrack
size_t ntrack
Number of along-track values.
Definition: libiasi.h:169
iasi_l2_t
IASI Level-2 data.
Definition: libiasi.h:198
iasi_l2_t::time
double time[L2_NTRACK][L2_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:204
iasi_l2_t::lon
double lon[L2_NTRACK][L2_NXTRACK]
Longitude [deg].
Definition: libiasi.h:210
iasi_l2_t::p
double p[L2_NLAY]
Pressure [hPa].
Definition: libiasi.h:216
iasi_l2_t::lat
double lat[L2_NTRACK][L2_NXTRACK]
Latitude [deg].
Definition: libiasi.h:213
iasi_l2_t::t
double t[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Temperature [K].
Definition: libiasi.h:219
iasi_l2_t::z
double z[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Geopotential height [km].
Definition: libiasi.h:207
iasi_l2_t::ntrack
size_t ntrack
Number of along-track values.
Definition: libiasi.h:201
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::Sat_lon
double Sat_lon[L1_NTRACK]
Estimated longitude of the satellite.
Definition: libiasi.h:312
iasi_rad_t::Sat_z
double Sat_z[L1_NTRACK]
Altitude of the satellite.
Definition: libiasi.h:309
iasi_rad_t::Rad
float Rad[L1_NTRACK][L1_NXTRACK][IASI_L1_NCHAN]
Radiance [W/(m^2 sr cm^-1)].
Definition: libiasi.h:306
iasi_rad_t::Sat_lat
double Sat_lat[L1_NTRACK]
Estimated latitude of the satellite.
Definition: libiasi.h:315
met_t
Meteorological data.
Definition: extract.c:86
met_t::h2o
float h2o[EX][EY][EP]
Water vapor volume mixing ratio [1].
Definition: extract.c:134
met_t::w
float w[EX][EY][EP]
Vertical wind [hPa/s].
Definition: extract.c:128
met_t::nx
int nx
Number of longitudes.
Definition: extract.c:92
met_t::ny
int ny
Number of latitudes.
Definition: extract.c:95
met_t::pt
double pt[EX][EY]
Tropopause pressure [hPa].
Definition: extract.c:113
met_t::o3
float o3[EX][EY][EP]
Ozone volume mixing ratio [1].
Definition: extract.c:137
met_t::np
int np
Number of pressure levels.
Definition: extract.c:98
met_t::t
float t[EX][EY][EP]
Temperature [K].
Definition: extract.c:119
met_t::u
float u[EX][EY][EP]
Zonal wind [m/s].
Definition: extract.c:122
met_t::lon
double lon[EX]
Longitude [deg].
Definition: extract.c:101
met_t::z
float z[EX][EY][EP]
Geopotential height [km].
Definition: extract.c:116
met_t::v
float v[EX][EY][EP]
Meridional wind [m/s].
Definition: extract.c:125
met_t::pv
float pv[EX][EY][EP]
Potential vorticity [PVU].
Definition: extract.c:131
met_t::time
double time
Time [s].
Definition: extract.c:89
met_t::ps
double ps[EX][EY]
Surface pressure [hPa].
Definition: extract.c:110
met_t::lat
double lat[EY]
Latitude [deg].
Definition: extract.c:104
met_t::p
double p[EP]
Pressure [hPa].
Definition: extract.c:107
Generated by doxygen 1.9.4