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