AIRS Code Collection
issifm.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 Dimensions...
30 ------------------------------------------------------------ */
31
33#define NZ 150
34
36#define NLON 4500
37
39#define NLAT 2400
40
41/* ------------------------------------------------------------
42 Structs...
43 ------------------------------------------------------------ */
44
46typedef struct {
47
49 int nz;
50
52 int nlon;
53
55 int nlat;
56
58 double lon[NLON];
59
61 double lat[NLAT];
62
64 float ps[NLON][NLAT];
65
67 float p[NLON][NLAT][NZ];
68
70 float t[NLON][NLAT][NZ];
71
73 float z[NLON][NLAT][NZ];
74
75} model_t;
76
77/* ------------------------------------------------------------
78 Functions...
79 ------------------------------------------------------------ */
80
82void intpol(
83 model_t * model,
84 double z,
85 double lon,
86 double lat,
87 double *p,
88 double *t);
89
91void smooth(
92 model_t * model);
93
95void write_nc(
96 char *filename,
97 wave_t * wave);
98
99/* ------------------------------------------------------------
100 Main...
101 ------------------------------------------------------------ */
102
103int main(
104 int argc,
105 char *argv[]) {
106
107 static ctl_t ctl;
108
109 static char kernel[LEN], pertname[LEN];
110
111 const double var_dh = 100.;
112
113 static double hyam[NZ], hybm[NZ], kz[NSHAPE], kw[NSHAPE];
114
115 static float *help;
116
117 static int init, ncid, dimid, varid, track0, track1, nk;
118
119 static size_t rs;
120
121 atm_t *atm;
122
123 obs_t *obs;
124
125 pert_t *pert;
126
127 wave_t *wave;
128
129 model_t *model;
130
131 /* ------------------------------------------------------------
132 Get control parameters...
133 ------------------------------------------------------------ */
134
135 /* Check arguments... */
136 if (argc < 8)
137 ERRMSG("Give parameters: <ctl> <model> <model.nc> <pert.nc>"
138 " <wave_airs.tab> <wave_model.tab> <wave_airs.nc> <wave_model.nc>");
139
140 /* Read control parameters... */
141 read_ctl(argc, argv, &ctl);
142 scan_ctl(argc, argv, "PERTNAME", -1, "4mu", pertname);
143 scan_ctl(argc, argv, "KERNEL", -1, "-", kernel);
144 const int extpol = (int) scan_ctl(argc, argv, "EXTPOL", -1, "0", NULL);
145 const int slant = (int) scan_ctl(argc, argv, "SLANT", -1, "1", NULL);
146 const double t_ovp = scan_ctl(argc, argv, "T_OVP", -1, "", NULL);
147
148 /* Set control parameters... */
149 ctl.write_bbt = 1;
150
151 /* Initialize look-up tables... */
152 tbl_t *tbl = read_tbl(&ctl);
153
154 /* ------------------------------------------------------------
155 Read model data...
156 ------------------------------------------------------------ */
157
158 /* Allocate... */
159 ALLOC(model, model_t, 1);
160 ALLOC(help, float,
161 NLON * NLAT * NZ);
162
163 /* Open file... */
164 printf("Read %s data: %s\n", argv[2], argv[3]);
165 NC(nc_open(argv[3], NC_NOWRITE, &ncid));
166
167 /* Check time... */
168 LOG(2, "Check time...");
169 if (nc_inq_dimid(ncid, "time", &dimid) != NC_NOERR)
170 NC(nc_inq_dimid(ncid, "time", &dimid));
171 NC(nc_inq_dimlen(ncid, dimid, &rs));
172 if (rs != 1)
173 ERRMSG("Only one time step is allowed!");
174
175 /* Read latitudes... */
176 LOG(2, "Read latitudes...");
177 if (nc_inq_dimid(ncid, "lat", &dimid) != NC_NOERR)
178 NC(nc_inq_dimid(ncid, "latitude", &dimid));
179 NC(nc_inq_dimlen(ncid, dimid, &rs));
180 model->nlat = (int) rs;
181 if (model->nlat > NLAT)
182 ERRMSG("Too many latitudes!");
183 if (nc_inq_varid(ncid, "lat", &varid) != NC_NOERR)
184 NC(nc_inq_varid(ncid, "latitude", &varid));
185 NC(nc_get_var_double(ncid, varid, model->lat));
186
187 /* Read longitudes... */
188 LOG(2, "Read longitudes...");
189 if (nc_inq_dimid(ncid, "lon", &dimid) != NC_NOERR)
190 NC(nc_inq_dimid(ncid, "longitude", &dimid));
191 NC(nc_inq_dimlen(ncid, dimid, &rs));
192 model->nlon = (int) rs;
193 if (model->nlon > NLON)
194 ERRMSG("Too many longitudes!");
195 if (nc_inq_varid(ncid, "lon", &varid))
196 NC(nc_inq_varid(ncid, "longitude", &varid));
197 NC(nc_get_var_double(ncid, varid, model->lon));
198
199 /* Read ICON data... */
200 if (strcasecmp(argv[2], "icon") == 0) {
201
202 /* Get height levels... */
203 LOG(2, "Read levels...");
204 NC(nc_inq_dimid(ncid, "height", &dimid));
205 NC(nc_inq_dimlen(ncid, dimid, &rs));
206 model->nz = (int) rs;
207 if (model->nz > NZ)
208 ERRMSG("Too many altitudes!");
209
210 /* Read height... */
211 LOG(2, "Read height...");
212 NC(nc_inq_varid(ncid, "z_mc", &varid));
213 NC(nc_get_var_float(ncid, varid, help));
214 for (int ilon = 0; ilon < model->nlon; ilon++)
215 for (int ilat = 0; ilat < model->nlat; ilat++)
216 for (int iz = 0; iz < model->nz; iz++)
217 model->z[ilon][ilat][iz] =
218 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
219 1e3);
220
221 /* Read temperature... */
222 LOG(2, "Read temperature...");
223 NC(nc_inq_varid(ncid, "temp", &varid));
224 NC(nc_get_var_float(ncid, varid, help));
225 for (int ilon = 0; ilon < model->nlon; ilon++)
226 for (int ilat = 0; ilat < model->nlat; ilat++)
227 for (int iz = 0; iz < model->nz; iz++)
228 model->t[ilon][ilat][iz] =
229 help[(iz * model->nlat + ilat) * model->nlon + ilon];
230
231 /* Read pressure... */
232 LOG(2, "Read pressure...");
233 NC(nc_inq_varid(ncid, "pres", &varid));
234 NC(nc_get_var_float(ncid, varid, help));
235 for (int ilon = 0; ilon < model->nlon; ilon++)
236 for (int ilat = 0; ilat < model->nlat; ilat++)
237 for (int iz = 0; iz < model->nz; iz++)
238 model->p[ilon][ilat][iz] =
239 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
240 1e2);
241 }
242
243 /* Read IFS data... */
244 else if (strcasecmp(argv[2], "ifs") == 0) {
245
246 /* Get height levels... */
247 LOG(2, "Read levels...");
248 NC(nc_inq_dimid(ncid, "lev_2", &dimid));
249 NC(nc_inq_dimlen(ncid, dimid, &rs));
250 model->nz = (int) rs;
251 if (model->nz > NZ)
252 ERRMSG("Too many altitudes!");
253
254 /* Read height... */
255 LOG(2, "Read height...");
256 NC(nc_inq_varid(ncid, "gh", &varid));
257 NC(nc_get_var_float(ncid, varid, help));
258 for (int ilon = 0; ilon < model->nlon; ilon++)
259 for (int ilat = 0; ilat < model->nlat; ilat++)
260 for (int iz = 0; iz < model->nz; iz++)
261 model->z[ilon][ilat][iz] =
262 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
263 1e3);
264
265 /* Read temperature... */
266 LOG(2, "Read temperature...");
267 NC(nc_inq_varid(ncid, "t", &varid));
268 NC(nc_get_var_float(ncid, varid, help));
269 for (int ilon = 0; ilon < model->nlon; ilon++)
270 for (int ilat = 0; ilat < model->nlat; ilat++)
271 for (int iz = 0; iz < model->nz; iz++)
272 model->t[ilon][ilat][iz] =
273 help[(iz * model->nlat + ilat) * model->nlon + ilon];
274
275 /* Read surface pressure... */
276 LOG(2, "Read surface pressure...");
277 NC(nc_inq_varid(ncid, "lnsp", &varid));
278 NC(nc_get_var_float(ncid, varid, help));
279 for (int ilon = 0; ilon < model->nlon; ilon++)
280 for (int ilat = 0; ilat < model->nlat; ilat++)
281 model->ps[ilon][ilat] = (float) exp(help[ilat * model->nlon + ilon]);
282
283 /* Read grid coefficients... */
284 LOG(2, "Read grid coefficients...");
285 NC(nc_inq_varid(ncid, "hyam", &varid));
286 NC(nc_get_var_double(ncid, varid, hyam));
287 NC(nc_inq_varid(ncid, "hybm", &varid));
288 NC(nc_get_var_double(ncid, varid, hybm));
289
290 /* Calculate pressure... */
291 LOG(2, "Calculate pressure...");
292 for (int ilon = 0; ilon < model->nlon; ilon++)
293 for (int ilat = 0; ilat < model->nlat; ilat++)
294 for (int iz = 0; iz < model->nz; iz++)
295 model->p[ilon][ilat][iz]
296 = (float) ((hyam[iz] + hybm[iz] * model->ps[ilon][ilat]) / 100.);
297 }
298
299 /* Read UM data... */
300 else if (strcasecmp(argv[2], "um") == 0) {
301
302 /* Get height levels... */
303 LOG(2, "Read levels...");
304 if (nc_inq_dimid(ncid, "RHO_TOP_eta_rho", &dimid) != NC_NOERR)
305 NC(nc_inq_dimid(ncid, "RHO_eta_rho", &dimid));
306 NC(nc_inq_dimlen(ncid, dimid, &rs));
307 model->nz = (int) rs;
308 if (model->nz > NZ)
309 ERRMSG("Too many altitudes!");
310
311 /* Read height... */
312 LOG(2, "Read height...");
313 if (nc_inq_varid(ncid, "STASH_m01s15i102_2", &varid) != NC_NOERR)
314 NC(nc_inq_varid(ncid, "STASH_m01s15i102", &varid));
315 NC(nc_get_var_float(ncid, varid, help));
316 for (int ilon = 0; ilon < model->nlon; ilon++)
317 for (int ilat = 0; ilat < model->nlat; ilat++)
318 for (int iz = 0; iz < model->nz; iz++)
319 model->z[ilon][ilat][iz] =
320 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
321 1e3);
322
323 /* Read temperature... */
324 LOG(2, "Read temperature...");
325 NC(nc_inq_varid(ncid, "STASH_m01s30i004", &varid));
326 NC(nc_get_var_float(ncid, varid, help));
327 for (int ilon = 0; ilon < model->nlon; ilon++)
328 for (int ilat = 0; ilat < model->nlat; ilat++)
329 for (int iz = 0; iz < model->nz; iz++)
330 model->t[ilon][ilat][iz] =
331 help[(iz * model->nlat + ilat) * model->nlon + ilon];
332
333 /* Read pressure... */
334 LOG(2, "Read pressure...");
335 NC(nc_inq_varid(ncid, "STASH_m01s00i407", &varid));
336 NC(nc_get_var_float(ncid, varid, help));
337 for (int ilon = 0; ilon < model->nlon; ilon++)
338 for (int ilat = 0; ilat < model->nlat; ilat++)
339 for (int iz = 0; iz < model->nz; iz++)
340 model->p[ilon][ilat][iz] =
341 0.01f * help[(iz * model->nlat + ilat) * model->nlon + ilon];
342 }
343
344 /* Read WRF data... */
345 else if (strcasecmp(argv[2], "wrf") == 0) {
346
347 /* Get height levels... */
348 LOG(2, "Read levels...");
349 NC(nc_inq_dimid(ncid, "bottom_top", &dimid));
350 NC(nc_inq_dimlen(ncid, dimid, &rs));
351 model->nz = (int) rs;
352 if (model->nz > NZ)
353 ERRMSG("Too many altitudes!");
354
355 /* Read height... */
356 LOG(2, "Read height...");
357 NC(nc_inq_varid(ncid, "z", &varid));
358 NC(nc_get_var_float(ncid, varid, help));
359 for (int ilon = 0; ilon < model->nlon; ilon++)
360 for (int ilat = 0; ilat < model->nlat; ilat++)
361 for (int iz = 0; iz < model->nz; iz++)
362 model->z[ilon][ilat][iz] =
363 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
364 1e3);
365
366 /* Read temperature... */
367 LOG(2, "Read temperature...");
368 NC(nc_inq_varid(ncid, "tk", &varid));
369 NC(nc_get_var_float(ncid, varid, help));
370 for (int ilon = 0; ilon < model->nlon; ilon++)
371 for (int ilat = 0; ilat < model->nlat; ilat++)
372 for (int iz = 0; iz < model->nz; iz++)
373 model->t[ilon][ilat][iz] =
374 help[(iz * model->nlat + ilat) * model->nlon + ilon];
375
376 /* Read pressure... */
377 LOG(2, "Read pressure...");
378 NC(nc_inq_varid(ncid, "p", &varid));
379 NC(nc_get_var_float(ncid, varid, help));
380 for (int ilon = 0; ilon < model->nlon; ilon++)
381 for (int ilat = 0; ilat < model->nlat; ilat++)
382 for (int iz = 0; iz < model->nz; iz++)
383 model->p[ilon][ilat][iz] =
384 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
385 1e2);
386 }
387
388 else
389 ERRMSG("Model type not supported!");
390
391 /* Close file... */
392 NC(nc_close(ncid));
393
394 /* Free... */
395 free(help);
396
397 /* Check data... */
398 LOG(2, "Checking...");
399 for (int ilon = 0; ilon < model->nlon; ilon++)
400 for (int ilat = 0; ilat < model->nlat; ilat++)
401 for (int iz = 0; iz < model->nz; iz++)
402 if (model->t[ilon][ilat][iz] <= 100
403 || model->t[ilon][ilat][iz] >= 400) {
404 model->p[ilon][ilat][iz] = GSL_NAN;
405 model->t[ilon][ilat][iz] = GSL_NAN;
406 model->z[ilon][ilat][iz] = GSL_NAN;
407 }
408
409 /* Smoothing of model data... */
410 LOG(2, "Smoothing...");
411 smooth(model);
412
413 /* Write info... */
414 for (int iz = 0; iz < model->nz; iz++)
415 printf("section_height: %d %g %g %g %g %g\n", iz,
416 model->z[model->nlon / 2][model->nlat / 2][iz],
417 model->lon[model->nlon / 2], model->lat[model->nlat / 2],
418 model->p[model->nlon / 2][model->nlat / 2][iz],
419 model->t[model->nlon / 2][model->nlat / 2][iz]);
420 for (int ilon = 0; ilon < model->nlon; ilon++)
421 printf("section_west_east: %d %g %g %g %g %g\n", ilon,
422 model->z[ilon][model->nlat / 2][model->nz / 2],
423 model->lon[ilon], model->lat[model->nlat / 2],
424 model->p[ilon][model->nlat / 2][model->nz / 2],
425 model->t[ilon][model->nlat / 2][model->nz / 2]);
426 for (int ilat = 0; ilat < model->nlat; ilat++)
427 printf("section_north_south: %d %g %g %g %g %g\n", ilat,
428 model->z[model->nlon / 2][ilat][model->nz / 2],
429 model->lon[model->nlon / 2], model->lat[ilat],
430 model->p[model->nlon / 2][ilat][model->nz / 2],
431 model->t[model->nlon / 2][ilat][model->nz / 2]);
432
433 /* ------------------------------------------------------------
434 Read AIRS perturbation data...
435 ------------------------------------------------------------ */
436
437 /* Allocate... */
438 ALLOC(atm, atm_t, 1);
439 ALLOC(obs, obs_t, 1);
440 ALLOC(pert, pert_t, 1);
441 ALLOC(wave, wave_t, 1);
442
443 /* Read perturbation data... */
444 read_pert(argv[4], pertname, pert);
445
446 /* Find track range... */
447 for (int itrack = 0; itrack < pert->ntrack; itrack++) {
448 if (pert->time[itrack][44] < t_ovp - 720 || itrack == 0)
449 track0 = itrack;
450 track1 = itrack;
451 if (pert->time[itrack][44] > t_ovp + 720)
452 break;
453 }
454
455 /* Convert to wave analysis struct... */
456 pert2wave(pert, wave, track0, track1, 0, pert->nxtrack - 1);
457
458 /* Estimate background... */
459 background_poly(wave, 5, 0);
460
461 /* Compute variance... */
462 variance(wave, var_dh);
463
464 /* Write observation wave struct... */
465 write_wave(argv[5], wave);
466 write_nc(argv[7], wave);
467
468 /* ------------------------------------------------------------
469 Run forward model...
470 ------------------------------------------------------------ */
471
472 /* Loop over AIRS geolocations... */
473 for (int itrack = track0; itrack <= track1; itrack++)
474 for (int ixtrack = 0; ixtrack < pert->nxtrack; ixtrack++) {
475
476 /* Write info... */
477 if (ixtrack == 0)
478 printf("Compute track %d / %d ...\n", itrack - track0 + 1,
479 track1 - track0 + 1);
480
481 /* Set observation data... */
482 obs->nr = 1;
483 obs->obsz[0] = 705;
484 obs->obslon[0] = pert->lon[itrack][44];
485 obs->obslat[0] = pert->lat[itrack][44];
486 obs->vpz[0] = 0;
487 obs->vplon[0] = pert->lon[itrack][ixtrack];
488 obs->vplat[0] = pert->lat[itrack][ixtrack];
489
490 /* Get Cartesian coordinates... */
491 double xm[3], xo[3], xs[3];
492 geo2cart(obs->obsz[0], obs->obslon[0], obs->obslat[0], xo);
493 geo2cart(obs->vpz[0], obs->vplon[0], obs->vplat[0], xs);
494
495 /* Set profile for atmospheric data... */
496 if (slant) {
497 atm->np = 0;
498 for (double f = 0.0; f <= 1.0; f += 0.0002) {
499 xm[0] = f * xo[0] + (1 - f) * xs[0];
500 xm[1] = f * xo[1] + (1 - f) * xs[1];
501 xm[2] = f * xo[2] + (1 - f) * xs[2];
502 cart2geo(xm, &atm->z[atm->np], &atm->lon[atm->np],
503 &atm->lat[atm->np]);
504 atm->time[atm->np] = pert->time[itrack][ixtrack];
505 if (atm->z[atm->np] < 10)
506 continue;
507 else if (atm->z[atm->np] > 90)
508 break;
509 else if ((++atm->np) >= NP)
510 ERRMSG("Too many ray path data points!");
511 }
512 } else {
513 atm->np = 0;
514 for (double f = 10.0; f <= 90.0; f += 0.2) {
515 atm->time[atm->np] = pert->time[itrack][ixtrack];
516 atm->z[atm->np] = f;
517 atm->lon[atm->np] = pert->lon[itrack][ixtrack];
518 atm->lat[atm->np] = pert->lat[itrack][ixtrack];
519 if ((++atm->np) >= NP)
520 ERRMSG("Too many ray path data points!");
521 }
522 }
523
524 /* Initialize with climatological data... */
525 climatology(&ctl, atm);
526
527 /* Interpolate model data... */
528 for (int ip = 0; ip < atm->np; ip++)
529 intpol(model, atm->z[ip], atm->lon[ip], atm->lat[ip],
530 &atm->p[ip], &atm->t[ip]);
531
532 /* Check profile... */
533 int okay = 1;
534 for (int ip = 0; ip < atm->np; ip++)
535 if (!gsl_finite(atm->p[ip]) || !gsl_finite(atm->t[ip]))
536 okay = 0;
537 if (!okay)
538 pert->bt[itrack][ixtrack] = GSL_NAN;
539 else {
540
541 /* Use kernel function... */
542 if (kernel[0] != '-') {
543
544 /* Read kernel function... */
545 if (!init) {
546 init = 1;
547 read_shape(kernel, kz, kw, &nk);
548 if (kz[0] > kz[1])
549 ERRMSG("Kernel function must be ascending!");
550 }
551
552 /* Calculate mean temperature... */
553 double wsum = 0;
554 pert->bt[itrack][ixtrack] = 0;
555 for (int ip = 0; ip < atm->np; ip++)
556 if (atm->z[ip] >= kz[0] && atm->z[ip] <= kz[nk - 1]) {
557 const int iz = locate_irr(kz, nk, atm->z[ip]);
558 const double w =
559 LIN(kz[iz], kw[iz], kz[iz + 1], kw[iz + 1], atm->z[ip]);
560 pert->bt[itrack][ixtrack] += w * atm->t[ip];
561 wsum += w;
562 }
563 pert->bt[itrack][ixtrack] /= wsum;
564 }
565
566 /* Use radiative transfer model... */
567 else {
568
569 /* Run forward model... */
570 formod(&ctl, tbl, atm, obs);
571
572 /* Get mean brightness temperature... */
573 pert->bt[itrack][ixtrack] = 0;
574 for (int id = 0; id < ctl.nd; id++)
575 pert->bt[itrack][ixtrack] += obs->rad[id][0] / ctl.nd;
576 }
577 }
578 }
579
580 /* Extrapolate... */
581 if (extpol)
582 for (int itrack = track0; itrack <= track1; itrack++) {
583 for (int ixtrack = 1; ixtrack < pert->nxtrack; ixtrack++)
584 if (!gsl_finite(pert->bt[itrack][ixtrack]))
585 pert->bt[itrack][ixtrack] = pert->bt[itrack][ixtrack - 1];
586 for (int ixtrack = pert->nxtrack - 2; ixtrack >= 0; ixtrack--)
587 if (!gsl_finite(pert->bt[itrack][ixtrack]))
588 pert->bt[itrack][ixtrack] = pert->bt[itrack][ixtrack + 1];
589 }
590
591 /* ------------------------------------------------------------
592 Write model perturbations...
593 ------------------------------------------------------------ */
594
595 /* Convert to wave analysis struct... */
596 pert2wave(pert, wave, track0, track1, 0, pert->nxtrack - 1);
597
598 /* Estimate background... */
599 background_poly(wave, 5, 0);
600
601 /* Compute variance... */
602 variance(wave, var_dh);
603
604 /* Write observation wave struct... */
605 write_wave(argv[6], wave);
606 write_nc(argv[8], wave);
607
608 /* Free... */
609 free(atm);
610 free(obs);
611 free(pert);
612 free(wave);
613 free(tbl);
614
615 return EXIT_SUCCESS;
616}
617
618/************************************************************************/
619
620void intpol(
621 model_t *model,
622 double z,
623 double lon,
624 double lat,
625 double *p,
626 double *t) {
627
628 double zd[NZ];
629
630 int iz;
631
632 /* Adjust longitude... */
633 if (model->lon[model->nlon - 1] > 180)
634 if (lon < 0)
635 lon += 360;
636
637 /* Check horizontal range... */
638 if (lon < model->lon[0]
639 || lon > model->lon[model->nlon - 1]
640 || lat < GSL_MIN(model->lat[0], model->lat[model->nlat - 1])
641 || lat > GSL_MAX(model->lat[0], model->lat[model->nlat - 1])) {
642 *p = GSL_NAN;
643 *t = GSL_NAN;
644 return;
645 }
646
647 /* Get indices... */
648 int ilon = locate_irr(model->lon, model->nlon, lon);
649 int ilat = locate_irr(model->lat, model->nlat, lat);
650
651 /* Check data... */
652 if (!gsl_finite(model->z[ilon][ilat][0])
653 || !gsl_finite(model->z[ilon][ilat][model->nz - 1])
654 || !gsl_finite(model->z[ilon][ilat + 1][model->nz - 1])
655 || !gsl_finite(model->z[ilon][ilat + 1][model->nz - 1])
656 || !gsl_finite(model->z[ilon + 1][ilat][model->nz - 1])
657 || !gsl_finite(model->z[ilon + 1][ilat][model->nz - 1])
658 || !gsl_finite(model->z[ilon + 1][ilat + 1][model->nz - 1])
659 || !gsl_finite(model->z[ilon + 1][ilat + 1][model->nz - 1])) {
660 *p = GSL_NAN;
661 *t = GSL_NAN;
662 return;
663 }
664
665 /* Check vertical range... */
666 if (z >
667 GSL_MAX(model->z[ilon][ilat][0], model->z[ilon][ilat][model->nz - 1])
668 || z < GSL_MIN(model->z[ilon][ilat][0],
669 model->z[ilon][ilat][model->nz - 1])
670 || z > GSL_MAX(model->z[ilon][ilat + 1][0],
671 model->z[ilon][ilat + 1][model->nz - 1])
672 || z < GSL_MIN(model->z[ilon][ilat + 1][0],
673 model->z[ilon][ilat + 1][model->nz - 1])
674 || z > GSL_MAX(model->z[ilon + 1][ilat][0],
675 model->z[ilon + 1][ilat][model->nz - 1])
676 || z < GSL_MIN(model->z[ilon + 1][ilat][0],
677 model->z[ilon + 1][ilat][model->nz - 1])
678 || z > GSL_MAX(model->z[ilon + 1][ilat + 1][0],
679 model->z[ilon + 1][ilat + 1][model->nz - 1])
680 || z < GSL_MIN(model->z[ilon + 1][ilat + 1][0],
681 model->z[ilon + 1][ilat + 1][model->nz - 1]))
682 return;
683
684 /* Interpolate vertically... */
685 for (iz = 0; iz < model->nz; iz++)
686 zd[iz] = model->z[ilon][ilat][iz];
687 iz = locate_irr(zd, model->nz, z);
688 double p00 = LIN(model->z[ilon][ilat][iz], model->p[ilon][ilat][iz],
689 model->z[ilon][ilat][iz + 1], model->p[ilon][ilat][iz + 1],
690 z);
691 double t00 = LIN(model->z[ilon][ilat][iz], model->t[ilon][ilat][iz],
692 model->z[ilon][ilat][iz + 1], model->t[ilon][ilat][iz + 1],
693 z);
694
695 for (iz = 0; iz < model->nz; iz++)
696 zd[iz] = model->z[ilon][ilat + 1][iz];
697 iz = locate_irr(zd, model->nz, z);
698 double p01 = LIN(model->z[ilon][ilat + 1][iz], model->p[ilon][ilat + 1][iz],
699 model->z[ilon][ilat + 1][iz + 1],
700 model->p[ilon][ilat + 1][iz + 1], z);
701 double t01 = LIN(model->z[ilon][ilat + 1][iz], model->t[ilon][ilat + 1][iz],
702 model->z[ilon][ilat + 1][iz + 1],
703 model->t[ilon][ilat + 1][iz + 1],
704 z);
705
706 for (iz = 0; iz < model->nz; iz++)
707 zd[iz] = model->z[ilon + 1][ilat][iz];
708 iz = locate_irr(zd, model->nz, z);
709 double p10 = LIN(model->z[ilon + 1][ilat][iz], model->p[ilon + 1][ilat][iz],
710 model->z[ilon + 1][ilat][iz + 1],
711 model->p[ilon + 1][ilat][iz + 1], z);
712 double t10 = LIN(model->z[ilon + 1][ilat][iz], model->t[ilon + 1][ilat][iz],
713 model->z[ilon + 1][ilat][iz + 1],
714 model->t[ilon + 1][ilat][iz + 1],
715 z);
716
717 for (iz = 0; iz < model->nz; iz++)
718 zd[iz] = model->z[ilon + 1][ilat + 1][iz];
719 iz = locate_irr(zd, model->nz, z);
720 double p11 =
721 LIN(model->z[ilon + 1][ilat + 1][iz], model->p[ilon + 1][ilat + 1][iz],
722 model->z[ilon + 1][ilat + 1][iz + 1],
723 model->p[ilon + 1][ilat + 1][iz + 1], z);
724 double t11 =
725 LIN(model->z[ilon + 1][ilat + 1][iz], model->t[ilon + 1][ilat + 1][iz],
726 model->z[ilon + 1][ilat + 1][iz + 1],
727 model->t[ilon + 1][ilat + 1][iz + 1], z);
728
729 /* Interpolate horizontally... */
730 p00 = LIN(model->lon[ilon], p00, model->lon[ilon + 1], p10, lon);
731 p11 = LIN(model->lon[ilon], p01, model->lon[ilon + 1], p11, lon);
732 *p = LIN(model->lat[ilat], p00, model->lat[ilat + 1], p11, lat);
733
734 t00 = LIN(model->lon[ilon], t00, model->lon[ilon + 1], t10, lon);
735 t11 = LIN(model->lon[ilon], t01, model->lon[ilon + 1], t11, lon);
736 *t = LIN(model->lat[ilat], t00, model->lat[ilat + 1], t11, lat);
737}
738
739/************************************************************************/
740
741void smooth(
742 model_t *model) {
743
744#define DLAT 14
745#define DLON 14
746
747 static float hp[NLON][NLAT], ht[NLON][NLAT], hz[NLON][NLAT];
748
749 const double fwhm = 20.0;
750
751 static double wy[DLAT + 1];
752
753 /* Set weights... */
754 const double dy = RE * M_PI / 180. * fabs(model->lat[1] - model->lat[0]);
755 for (int ilat = 0; ilat <= DLAT; ilat++)
756 wy[ilat] = exp(-0.5 * POW2(ilat * dy * 2.35482 / fwhm));
757
758 /* Loop over height levels... */
759 for (int iz = 0; iz < model->nz; iz++) {
760
761 /* Write info... */
762 printf("Smoothing level %d / %d ...\n", iz + 1, model->nz);
763
764 /* Copy data... */
765#pragma omp parallel for collapse(2)
766 for (int ilon = 0; ilon < model->nlon; ilon++)
767 for (int ilat = 0; ilat < model->nlat; ilat++) {
768 hp[ilon][ilat] = model->p[ilon][ilat][iz];
769 ht[ilon][ilat] = model->t[ilon][ilat][iz];
770 hz[ilon][ilat] = model->z[ilon][ilat][iz];
771 }
772
773 /* Loop over latitudes... */
774#pragma omp parallel
775 {
776 /* Thread-private wx */
777 double wx[DLON + 1];
778
779#pragma omp for
780 for (int ilat = 0; ilat < model->nlat; ilat++) {
781
782 /* Set weights... */
783 const double dx =
784 RE * M_PI / 180. * cos(model->lat[ilat] * M_PI / 180.) *
785 fabs(model->lon[1] - model->lon[0]);
786
787 for (int ilon = 0; ilon <= DLON; ilon++)
788 wx[ilon] = exp(-0.5 * POW2(ilon * dx * 2.35482 / fwhm));
789
790 /* Loop over longitudes... */
791 for (int ilon = 0; ilon < model->nlon; ilon++) {
792 float wsum = 0;
793 model->p[ilon][ilat][iz] = 0;
794 model->t[ilon][ilat][iz] = 0;
795 model->z[ilon][ilat][iz] = 0;
796
797 const int ilo0 = (ilon - DLON > 0) ? ilon - DLON : 0;
798 const int ilo1 =
799 (ilon + DLON < model->nlon - 1) ? ilon + DLON : model->nlon - 1;
800 const int ila0 = (ilat - DLAT > 0) ? ilat - DLAT : 0;
801 const int ila1 =
802 (ilat + DLAT < model->nlat - 1) ? ilat + DLAT : model->nlat - 1;
803
804 for (int ilon2 = ilo0; ilon2 <= ilo1; ++ilon2) {
805 const double wxh = wx[abs(ilon2 - ilon)];
806 for (int ilat2 = ila0; ilat2 <= ila1; ++ilat2) {
807 const float w = (float) (wxh * wy[abs(ilat2 - ilat)]);
808 model->p[ilon][ilat][iz] += w * hp[ilon2][ilat2];
809 model->t[ilon][ilat][iz] += w * ht[ilon2][ilat2];
810 model->z[ilon][ilat][iz] += w * hz[ilon2][ilat2];
811 wsum += w;
812 }
813 }
814
815 model->p[ilon][ilat][iz] /= wsum;
816 model->t[ilon][ilat][iz] /= wsum;
817 model->z[ilon][ilat][iz] /= wsum;
818 }
819 }
820 }
821 }
822}
823
824/************************************************************************/
825
826void write_nc(
827 char *filename,
828 wave_t *wave) {
829
830 static double help[WX * WY];
831
832 int ncid, dimid[10], lon_id, lat_id, bt_id, pt_id, var_id;
833
834 /* Create netCDF file... */
835 NC(nc_create(filename, NC_CLOBBER, &ncid));
836
837 /* Set dimensions... */
838 NC(nc_def_dim(ncid, "NTRACK", (size_t) wave->ny, &dimid[0]));
839 NC(nc_def_dim(ncid, "NXTRACK", (size_t) wave->nx, &dimid[1]));
840
841 /* Add variables... */
842 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 2, dimid, &lon_id));
843 add_att(ncid, lon_id, "deg", "footprint longitude");
844 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 2, dimid, &lat_id));
845 add_att(ncid, lat_id, "deg", "footprint latitude");
846 NC(nc_def_var(ncid, "bt", NC_FLOAT, 2, dimid, &bt_id));
847 add_att(ncid, bt_id, "K", "brightness temperature");
848 NC(nc_def_var(ncid, "bt_pt", NC_FLOAT, 2, dimid, &pt_id));
849 add_att(ncid, pt_id, "K", "brightness temperature perturbation");
850 NC(nc_def_var(ncid, "bt_var", NC_FLOAT, 2, dimid, &var_id));
851 add_att(ncid, var_id, "K^2", "brightness temperature variance");
852
853 /* Leave define mode... */
854 NC(nc_enddef(ncid));
855
856 /* Write data... */
857 for (int ix = 0; ix < wave->nx; ix++)
858 for (int iy = 0; iy < wave->ny; iy++)
859 help[iy * wave->nx + ix] = wave->lon[ix][iy];
860 NC(nc_put_var_double(ncid, lon_id, help));
861 for (int ix = 0; ix < wave->nx; ix++)
862 for (int iy = 0; iy < wave->ny; iy++)
863 help[iy * wave->nx + ix] = wave->lat[ix][iy];
864 NC(nc_put_var_double(ncid, lat_id, help));
865 for (int ix = 0; ix < wave->nx; ix++)
866 for (int iy = 0; iy < wave->ny; iy++)
867 help[iy * wave->nx + ix] = wave->temp[ix][iy];
868 NC(nc_put_var_double(ncid, bt_id, help));
869 for (int ix = 0; ix < wave->nx; ix++)
870 for (int iy = 0; iy < wave->ny; iy++)
871 help[iy * wave->nx + ix] = wave->pt[ix][iy];
872 NC(nc_put_var_double(ncid, pt_id, help));
873 for (int ix = 0; ix < wave->nx; ix++)
874 for (int iy = 0; iy < wave->ny; iy++)
875 help[iy * wave->nx + ix] = wave->var[ix][iy];
876 NC(nc_put_var_double(ncid, var_id, help));
877
878 /* Close file... */
879 NC(nc_close(ncid));
880}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: diff_apr.c:35
main
int main(int argc, char *argv[])
Definition: issifm.c:103
intpol
void intpol(model_t *model, double z, double lon, double lat, double *p, double *t)
Interpolation of model data.
Definition: issifm.c:620
write_nc
void write_nc(char *filename, wave_t *wave)
Write wave struct to netCDF file.
Definition: issifm.c:826
DLAT
#define DLAT
NZ
#define NZ
Maximum number of model levels.
Definition: issifm.c:33
smooth
void smooth(model_t *model)
Smoothing of model data.
Definition: issifm.c:741
DLON
#define DLON
NLAT
#define NLAT
Maximum number of model latitudes.
Definition: issifm.c:39
NLON
#define NLON
Maximum number of model longitudes.
Definition: issifm.c:36
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:1534
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:949
read_pert
void read_pert(char *filename, char *pertname, pert_t *pert)
Read radiance perturbation data.
Definition: libairs.c:1087
write_wave
void write_wave(char *filename, wave_t *wave)
Write wave analysis data.
Definition: libairs.c:1691
libairs.h
AIRS Code Collection library declarations.
WX
#define WX
Across-track size of wave analysis data.
Definition: libairs.h:129
WY
#define WY
Along-track size of wave analysis data.
Definition: libairs.h:132
model_t
Model data.
Definition: issifm.c:46
model_t::lat
double lat[NLAT]
Latitude [deg].
Definition: issifm.c:61
model_t::t
float t[NLON][NLAT][NZ]
Temperature [K].
Definition: issifm.c:70
model_t::lon
double lon[NLON]
Longitude [deg].
Definition: issifm.c:58
model_t::ps
float ps[NLON][NLAT]
Surface pressure [hPa].
Definition: issifm.c:64
model_t::nz
int nz
Number of vertical levels.
Definition: issifm.c:49
model_t::z
float z[NLON][NLAT][NZ]
Height [km].
Definition: issifm.c:73
model_t::nlon
int nlon
Number of longitudes.
Definition: issifm.c:52
model_t::nlat
int nlat
Number of latitudes.
Definition: issifm.c:55
model_t::p
float p[NLON][NLAT][NZ]
Pressure [hPa].
Definition: issifm.c:67
pert_t
Perturbation data.
Definition: libairs.h:205
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libairs.h:214
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::temp
double temp[WX][WY]
Temperature [K].
Definition: libairs.h:314
wave_t::lon
double lon[WX][WY]
Longitude [deg].
Definition: libairs.h:302
wave_t::lat
double lat[WX][WY]
Latitude [deg].
Definition: libairs.h:305
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libairs.h:320
Generated by doxygen 1.9.4