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 /* ------------------------------------------------------------
153 Read model data...
154 ------------------------------------------------------------ */
155
156 /* Allocate... */
157 ALLOC(model, model_t, 1);
158 ALLOC(help, float,
159 NLON * NLAT * NZ);
160
161 /* Open file... */
162 printf("Read %s data: %s\n", argv[2], argv[3]);
163 NC(nc_open(argv[3], NC_NOWRITE, &ncid));
164
165 /* Check time... */
166 LOG(2, "Check time...");
167 if (nc_inq_dimid(ncid, "time", &dimid) != NC_NOERR)
168 NC(nc_inq_dimid(ncid, "time", &dimid));
169 NC(nc_inq_dimlen(ncid, dimid, &rs));
170 if (rs != 1)
171 ERRMSG("Only one time step is allowed!");
172
173 /* Read latitudes... */
174 LOG(2, "Read latitudes...");
175 if (nc_inq_dimid(ncid, "lat", &dimid) != NC_NOERR)
176 NC(nc_inq_dimid(ncid, "latitude", &dimid));
177 NC(nc_inq_dimlen(ncid, dimid, &rs));
178 model->nlat = (int) rs;
179 if (model->nlat > NLAT)
180 ERRMSG("Too many latitudes!");
181 if (nc_inq_varid(ncid, "lat", &varid) != NC_NOERR)
182 NC(nc_inq_varid(ncid, "latitude", &varid));
183 NC(nc_get_var_double(ncid, varid, model->lat));
184
185 /* Read longitudes... */
186 LOG(2, "Read longitudes...");
187 if (nc_inq_dimid(ncid, "lon", &dimid) != NC_NOERR)
188 NC(nc_inq_dimid(ncid, "longitude", &dimid));
189 NC(nc_inq_dimlen(ncid, dimid, &rs));
190 model->nlon = (int) rs;
191 if (model->nlon > NLON)
192 ERRMSG("Too many longitudes!");
193 if (nc_inq_varid(ncid, "lon", &varid))
194 NC(nc_inq_varid(ncid, "longitude", &varid));
195 NC(nc_get_var_double(ncid, varid, model->lon));
196
197 /* Read ICON data... */
198 if (strcasecmp(argv[2], "icon") == 0) {
199
200 /* Get height levels... */
201 LOG(2, "Read levels...");
202 NC(nc_inq_dimid(ncid, "height", &dimid));
203 NC(nc_inq_dimlen(ncid, dimid, &rs));
204 model->nz = (int) rs;
205 if (model->nz > NZ)
206 ERRMSG("Too many altitudes!");
207
208 /* Read height... */
209 LOG(2, "Read height...");
210 NC(nc_inq_varid(ncid, "z_mc", &varid));
211 NC(nc_get_var_float(ncid, varid, help));
212 for (int ilon = 0; ilon < model->nlon; ilon++)
213 for (int ilat = 0; ilat < model->nlat; ilat++)
214 for (int iz = 0; iz < model->nz; iz++)
215 model->z[ilon][ilat][iz] =
216 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
217 1e3);
218
219 /* Read temperature... */
220 LOG(2, "Read temperature...");
221 NC(nc_inq_varid(ncid, "temp", &varid));
222 NC(nc_get_var_float(ncid, varid, help));
223 for (int ilon = 0; ilon < model->nlon; ilon++)
224 for (int ilat = 0; ilat < model->nlat; ilat++)
225 for (int iz = 0; iz < model->nz; iz++)
226 model->t[ilon][ilat][iz] =
227 help[(iz * model->nlat + ilat) * model->nlon + ilon];
228
229 /* Read pressure... */
230 LOG(2, "Read pressure...");
231 NC(nc_inq_varid(ncid, "pres", &varid));
232 NC(nc_get_var_float(ncid, varid, help));
233 for (int ilon = 0; ilon < model->nlon; ilon++)
234 for (int ilat = 0; ilat < model->nlat; ilat++)
235 for (int iz = 0; iz < model->nz; iz++)
236 model->p[ilon][ilat][iz] =
237 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
238 1e2);
239 }
240
241 /* Read IFS data... */
242 else if (strcasecmp(argv[2], "ifs") == 0) {
243
244 /* Get height levels... */
245 LOG(2, "Read levels...");
246 NC(nc_inq_dimid(ncid, "lev_2", &dimid));
247 NC(nc_inq_dimlen(ncid, dimid, &rs));
248 model->nz = (int) rs;
249 if (model->nz > NZ)
250 ERRMSG("Too many altitudes!");
251
252 /* Read height... */
253 LOG(2, "Read height...");
254 NC(nc_inq_varid(ncid, "gh", &varid));
255 NC(nc_get_var_float(ncid, varid, help));
256 for (int ilon = 0; ilon < model->nlon; ilon++)
257 for (int ilat = 0; ilat < model->nlat; ilat++)
258 for (int iz = 0; iz < model->nz; iz++)
259 model->z[ilon][ilat][iz] =
260 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
261 1e3);
262
263 /* Read temperature... */
264 LOG(2, "Read temperature...");
265 NC(nc_inq_varid(ncid, "t", &varid));
266 NC(nc_get_var_float(ncid, varid, help));
267 for (int ilon = 0; ilon < model->nlon; ilon++)
268 for (int ilat = 0; ilat < model->nlat; ilat++)
269 for (int iz = 0; iz < model->nz; iz++)
270 model->t[ilon][ilat][iz] =
271 help[(iz * model->nlat + ilat) * model->nlon + ilon];
272
273 /* Read surface pressure... */
274 LOG(2, "Read surface pressure...");
275 NC(nc_inq_varid(ncid, "lnsp", &varid));
276 NC(nc_get_var_float(ncid, varid, help));
277 for (int ilon = 0; ilon < model->nlon; ilon++)
278 for (int ilat = 0; ilat < model->nlat; ilat++)
279 model->ps[ilon][ilat] = (float) exp(help[ilat * model->nlon + ilon]);
280
281 /* Read grid coefficients... */
282 LOG(2, "Read grid coefficients...");
283 NC(nc_inq_varid(ncid, "hyam", &varid));
284 NC(nc_get_var_double(ncid, varid, hyam));
285 NC(nc_inq_varid(ncid, "hybm", &varid));
286 NC(nc_get_var_double(ncid, varid, hybm));
287
288 /* Calculate pressure... */
289 LOG(2, "Calculate pressure...");
290 for (int ilon = 0; ilon < model->nlon; ilon++)
291 for (int ilat = 0; ilat < model->nlat; ilat++)
292 for (int iz = 0; iz < model->nz; iz++)
293 model->p[ilon][ilat][iz]
294 = (float) ((hyam[iz] + hybm[iz] * model->ps[ilon][ilat]) / 100.);
295 }
296
297 /* Read UM data... */
298 else if (strcasecmp(argv[2], "um") == 0) {
299
300 /* Get height levels... */
301 LOG(2, "Read levels...");
302 if (nc_inq_dimid(ncid, "RHO_TOP_eta_rho", &dimid) != NC_NOERR)
303 NC(nc_inq_dimid(ncid, "RHO_eta_rho", &dimid));
304 NC(nc_inq_dimlen(ncid, dimid, &rs));
305 model->nz = (int) rs;
306 if (model->nz > NZ)
307 ERRMSG("Too many altitudes!");
308
309 /* Read height... */
310 LOG(2, "Read height...");
311 if (nc_inq_varid(ncid, "STASH_m01s15i102_2", &varid) != NC_NOERR)
312 NC(nc_inq_varid(ncid, "STASH_m01s15i102", &varid));
313 NC(nc_get_var_float(ncid, varid, help));
314 for (int ilon = 0; ilon < model->nlon; ilon++)
315 for (int ilat = 0; ilat < model->nlat; ilat++)
316 for (int iz = 0; iz < model->nz; iz++)
317 model->z[ilon][ilat][iz] =
318 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
319 1e3);
320
321 /* Read temperature... */
322 LOG(2, "Read temperature...");
323 NC(nc_inq_varid(ncid, "STASH_m01s30i004", &varid));
324 NC(nc_get_var_float(ncid, varid, help));
325 for (int ilon = 0; ilon < model->nlon; ilon++)
326 for (int ilat = 0; ilat < model->nlat; ilat++)
327 for (int iz = 0; iz < model->nz; iz++)
328 model->t[ilon][ilat][iz] =
329 help[(iz * model->nlat + ilat) * model->nlon + ilon];
330
331 /* Read pressure... */
332 LOG(2, "Read pressure...");
333 NC(nc_inq_varid(ncid, "STASH_m01s00i407", &varid));
334 NC(nc_get_var_float(ncid, varid, help));
335 for (int ilon = 0; ilon < model->nlon; ilon++)
336 for (int ilat = 0; ilat < model->nlat; ilat++)
337 for (int iz = 0; iz < model->nz; iz++)
338 model->p[ilon][ilat][iz] =
339 0.01f * help[(iz * model->nlat + ilat) * model->nlon + ilon];
340 }
341
342 /* Read WRF data... */
343 else if (strcasecmp(argv[2], "wrf") == 0) {
344
345 /* Get height levels... */
346 LOG(2, "Read levels...");
347 NC(nc_inq_dimid(ncid, "bottom_top", &dimid));
348 NC(nc_inq_dimlen(ncid, dimid, &rs));
349 model->nz = (int) rs;
350 if (model->nz > NZ)
351 ERRMSG("Too many altitudes!");
352
353 /* Read height... */
354 LOG(2, "Read height...");
355 NC(nc_inq_varid(ncid, "z", &varid));
356 NC(nc_get_var_float(ncid, varid, help));
357 for (int ilon = 0; ilon < model->nlon; ilon++)
358 for (int ilat = 0; ilat < model->nlat; ilat++)
359 for (int iz = 0; iz < model->nz; iz++)
360 model->z[ilon][ilat][iz] =
361 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
362 1e3);
363
364 /* Read temperature... */
365 LOG(2, "Read temperature...");
366 NC(nc_inq_varid(ncid, "tk", &varid));
367 NC(nc_get_var_float(ncid, varid, help));
368 for (int ilon = 0; ilon < model->nlon; ilon++)
369 for (int ilat = 0; ilat < model->nlat; ilat++)
370 for (int iz = 0; iz < model->nz; iz++)
371 model->t[ilon][ilat][iz] =
372 help[(iz * model->nlat + ilat) * model->nlon + ilon];
373
374 /* Read pressure... */
375 LOG(2, "Read pressure...");
376 NC(nc_inq_varid(ncid, "p", &varid));
377 NC(nc_get_var_float(ncid, varid, help));
378 for (int ilon = 0; ilon < model->nlon; ilon++)
379 for (int ilat = 0; ilat < model->nlat; ilat++)
380 for (int iz = 0; iz < model->nz; iz++)
381 model->p[ilon][ilat][iz] =
382 (float) (help[(iz * model->nlat + ilat) * model->nlon + ilon] /
383 1e2);
384 }
385
386 else
387 ERRMSG("Model type not supported!");
388
389 /* Close file... */
390 NC(nc_close(ncid));
391
392 /* Free... */
393 free(help);
394
395 /* Check data... */
396 LOG(2, "Checking...");
397 for (int ilon = 0; ilon < model->nlon; ilon++)
398 for (int ilat = 0; ilat < model->nlat; ilat++)
399 for (int iz = 0; iz < model->nz; iz++)
400 if (model->t[ilon][ilat][iz] <= 100
401 || model->t[ilon][ilat][iz] >= 400) {
402 model->p[ilon][ilat][iz] = GSL_NAN;
403 model->t[ilon][ilat][iz] = GSL_NAN;
404 model->z[ilon][ilat][iz] = GSL_NAN;
405 }
406
407 /* Smoothing of model data... */
408 LOG(2, "Smoothing...");
409 smooth(model);
410
411 /* Write info... */
412 for (int iz = 0; iz < model->nz; iz++)
413 printf("section_height: %d %g %g %g %g %g\n", iz,
414 model->z[model->nlon / 2][model->nlat / 2][iz],
415 model->lon[model->nlon / 2], model->lat[model->nlat / 2],
416 model->p[model->nlon / 2][model->nlat / 2][iz],
417 model->t[model->nlon / 2][model->nlat / 2][iz]);
418 for (int ilon = 0; ilon < model->nlon; ilon++)
419 printf("section_west_east: %d %g %g %g %g %g\n", ilon,
420 model->z[ilon][model->nlat / 2][model->nz / 2],
421 model->lon[ilon], model->lat[model->nlat / 2],
422 model->p[ilon][model->nlat / 2][model->nz / 2],
423 model->t[ilon][model->nlat / 2][model->nz / 2]);
424 for (int ilat = 0; ilat < model->nlat; ilat++)
425 printf("section_north_south: %d %g %g %g %g %g\n", ilat,
426 model->z[model->nlon / 2][ilat][model->nz / 2],
427 model->lon[model->nlon / 2], model->lat[ilat],
428 model->p[model->nlon / 2][ilat][model->nz / 2],
429 model->t[model->nlon / 2][ilat][model->nz / 2]);
430
431 /* ------------------------------------------------------------
432 Read AIRS perturbation data...
433 ------------------------------------------------------------ */
434
435 /* Allocate... */
436 ALLOC(atm, atm_t, 1);
437 ALLOC(obs, obs_t, 1);
438 ALLOC(pert, pert_t, 1);
439 ALLOC(wave, wave_t, 1);
440
441 /* Read perturbation data... */
442 read_pert(argv[4], pertname, pert);
443
444 /* Find track range... */
445 for (int itrack = 0; itrack < pert->ntrack; itrack++) {
446 if (pert->time[itrack][44] < t_ovp - 720 || itrack == 0)
447 track0 = itrack;
448 track1 = itrack;
449 if (pert->time[itrack][44] > t_ovp + 720)
450 break;
451 }
452
453 /* Convert to wave analysis struct... */
454 pert2wave(pert, wave, track0, track1, 0, pert->nxtrack - 1);
455
456 /* Estimate background... */
457 background_poly(wave, 5, 0);
458
459 /* Compute variance... */
460 variance(wave, var_dh);
461
462 /* Write observation wave struct... */
463 write_wave(argv[5], wave);
464 write_nc(argv[7], wave);
465
466 /* ------------------------------------------------------------
467 Run forward model...
468 ------------------------------------------------------------ */
469
470 /* Loop over AIRS geolocations... */
471 for (int itrack = track0; itrack <= track1; itrack++)
472 for (int ixtrack = 0; ixtrack < pert->nxtrack; ixtrack++) {
473
474 /* Write info... */
475 if (ixtrack == 0)
476 printf("Compute track %d / %d ...\n", itrack - track0 + 1,
477 track1 - track0 + 1);
478
479 /* Set observation data... */
480 obs->nr = 1;
481 obs->obsz[0] = 705;
482 obs->obslon[0] = pert->lon[itrack][44];
483 obs->obslat[0] = pert->lat[itrack][44];
484 obs->vpz[0] = 0;
485 obs->vplon[0] = pert->lon[itrack][ixtrack];
486 obs->vplat[0] = pert->lat[itrack][ixtrack];
487
488 /* Get Cartesian coordinates... */
489 geo2cart(obs->obsz[0], obs->obslon[0], obs->obslat[0], xo);
490 geo2cart(obs->vpz[0], obs->vplon[0], obs->vplat[0], xs);
491
492 /* Set profile for atmospheric data... */
493 if (slant) {
494 atm->np = 0;
495 for (double f = 0.0; f <= 1.0; f += 0.0002) {
496 xm[0] = f * xo[0] + (1 - f) * xs[0];
497 xm[1] = f * xo[1] + (1 - f) * xs[1];
498 xm[2] = f * xo[2] + (1 - f) * xs[2];
499 cart2geo(xm, &atm->z[atm->np], &atm->lon[atm->np],
500 &atm->lat[atm->np]);
501 atm->time[atm->np] = pert->time[itrack][ixtrack];
502 if (atm->z[atm->np] < 10)
503 continue;
504 else if (atm->z[atm->np] > 90)
505 break;
506 else if ((++atm->np) >= NP)
507 ERRMSG("Too many altitudes!");
508 }
509 } else {
510 atm->np = 0;
511 for (double f = 10.0; f <= 90.0; f += 0.2) {
512 atm->time[atm->np] = pert->time[itrack][ixtrack];
513 atm->z[atm->np] = f;
514 atm->lon[atm->np] = pert->lon[itrack][ixtrack];
515 atm->lat[atm->np] = pert->lat[itrack][ixtrack];
516 if ((++atm->np) >= NP)
517 ERRMSG("Too many altitudes!");
518 }
519 }
520
521 /* Initialize with climatological data... */
522 climatology(&ctl, atm);
523
524 /* Interpolate model data... */
525 for (int ip = 0; ip < atm->np; ip++)
526 intpol(model, atm->z[ip], atm->lon[ip], atm->lat[ip],
527 &atm->p[ip], &atm->t[ip]);
528
529 /* Check profile... */
530 int okay = 1;
531 for (int ip = 0; ip < atm->np; ip++)
532 if (!gsl_finite(atm->p[ip]) || !gsl_finite(atm->t[ip]))
533 okay = 0;
534 if (!okay)
535 pert->bt[itrack][ixtrack] = GSL_NAN;
536 else {
537
538 /* Use kernel function... */
539 if (kernel[0] != '-') {
540
541 /* Read kernel function... */
542 if (!init) {
543 init = 1;
544 read_shape(kernel, kz, kw, &nk);
545 if (kz[0] > kz[1])
546 ERRMSG("Kernel function must be ascending!");
547 }
548
549 /* Calculate mean temperature... */
550 pert->bt[itrack][ixtrack] = wsum = 0;
551 for (int ip = 0; ip < atm->np; ip++)
552 if (atm->z[ip] >= kz[0] && atm->z[ip] <= kz[nk - 1]) {
553 const int iz = locate_irr(kz, nk, atm->z[ip]);
554 const double w =
555 LIN(kz[iz], kw[iz], kz[iz + 1], kw[iz + 1], atm->z[ip]);
556 pert->bt[itrack][ixtrack] += w * atm->t[ip];
557 wsum += w;
558 }
559 pert->bt[itrack][ixtrack] /= wsum;
560 }
561
562 /* Use radiative transfer model... */
563 else {
564
565 /* Run forward model... */
566 formod(&ctl, atm, obs);
567
568 /* Get mean brightness temperature... */
569 pert->bt[itrack][ixtrack] = 0;
570 for (int id = 0; id < ctl.nd; id++)
571 pert->bt[itrack][ixtrack] += obs->rad[id][0] / ctl.nd;
572 }
573 }
574 }
575
576 /* Extrapolate... */
577 if (extpol)
578 for (int itrack = track0; itrack <= track1; itrack++) {
579 for (int ixtrack = 1; ixtrack < pert->nxtrack; ixtrack++)
580 if (!gsl_finite(pert->bt[itrack][ixtrack]))
581 pert->bt[itrack][ixtrack] = pert->bt[itrack][ixtrack - 1];
582 for (int ixtrack = pert->nxtrack - 2; ixtrack >= 0; ixtrack--)
583 if (!gsl_finite(pert->bt[itrack][ixtrack]))
584 pert->bt[itrack][ixtrack] = pert->bt[itrack][ixtrack + 1];
585 }
586
587 /* ------------------------------------------------------------
588 Write model perturbations...
589 ------------------------------------------------------------ */
590
591 /* Convert to wave analysis struct... */
592 pert2wave(pert, wave, track0, track1, 0, pert->nxtrack - 1);
593
594 /* Estimate background... */
595 background_poly(wave, 5, 0);
596
597 /* Compute variance... */
598 variance(wave, var_dh);
599
600 /* Write observation wave struct... */
601 write_wave(argv[6], wave);
602 write_nc(argv[8], wave);
603
604 /* Free... */
605 free(atm);
606 free(obs);
607 free(pert);
608 free(wave);
609
610 return EXIT_SUCCESS;
611}
612
613/************************************************************************/
614
615void intpol(
616 model_t *model,
617 double z,
618 double lon,
619 double lat,
620 double *p,
621 double *t) {
622
623 double zd[NZ];
624
625 int iz;
626
627 /* Adjust longitude... */
628 if (model->lon[model->nlon - 1] > 180)
629 if (lon < 0)
630 lon += 360;
631
632 /* Check horizontal range... */
633 if (lon < model->lon[0]
634 || lon > model->lon[model->nlon - 1]
635 || lat < GSL_MIN(model->lat[0], model->lat[model->nlat - 1])
636 || lat > GSL_MAX(model->lat[0], model->lat[model->nlat - 1])) {
637 *p = GSL_NAN;
638 *t = GSL_NAN;
639 return;
640 }
641
642 /* Get indices... */
643 int ilon = locate_irr(model->lon, model->nlon, lon);
644 int ilat = locate_irr(model->lat, model->nlat, lat);
645
646 /* Check data... */
647 if (!gsl_finite(model->z[ilon][ilat][0])
648 || !gsl_finite(model->z[ilon][ilat][model->nz - 1])
649 || !gsl_finite(model->z[ilon][ilat + 1][model->nz - 1])
650 || !gsl_finite(model->z[ilon][ilat + 1][model->nz - 1])
651 || !gsl_finite(model->z[ilon + 1][ilat][model->nz - 1])
652 || !gsl_finite(model->z[ilon + 1][ilat][model->nz - 1])
653 || !gsl_finite(model->z[ilon + 1][ilat + 1][model->nz - 1])
654 || !gsl_finite(model->z[ilon + 1][ilat + 1][model->nz - 1])) {
655 *p = GSL_NAN;
656 *t = GSL_NAN;
657 return;
658 }
659
660 /* Check vertical range... */
661 if (z >
662 GSL_MAX(model->z[ilon][ilat][0], model->z[ilon][ilat][model->nz - 1])
663 || z < GSL_MIN(model->z[ilon][ilat][0],
664 model->z[ilon][ilat][model->nz - 1])
665 || z > GSL_MAX(model->z[ilon][ilat + 1][0],
666 model->z[ilon][ilat + 1][model->nz - 1])
667 || z < GSL_MIN(model->z[ilon][ilat + 1][0],
668 model->z[ilon][ilat + 1][model->nz - 1])
669 || z > GSL_MAX(model->z[ilon + 1][ilat][0],
670 model->z[ilon + 1][ilat][model->nz - 1])
671 || z < GSL_MIN(model->z[ilon + 1][ilat][0],
672 model->z[ilon + 1][ilat][model->nz - 1])
673 || z > GSL_MAX(model->z[ilon + 1][ilat + 1][0],
674 model->z[ilon + 1][ilat + 1][model->nz - 1])
675 || z < GSL_MIN(model->z[ilon + 1][ilat + 1][0],
676 model->z[ilon + 1][ilat + 1][model->nz - 1]))
677 return;
678
679 /* Interpolate vertically... */
680 for (iz = 0; iz < model->nz; iz++)
681 zd[iz] = model->z[ilon][ilat][iz];
682 iz = locate_irr(zd, model->nz, z);
683 double p00 = LIN(model->z[ilon][ilat][iz], model->p[ilon][ilat][iz],
684 model->z[ilon][ilat][iz + 1], model->p[ilon][ilat][iz + 1],
685 z);
686 double t00 = LIN(model->z[ilon][ilat][iz], model->t[ilon][ilat][iz],
687 model->z[ilon][ilat][iz + 1], model->t[ilon][ilat][iz + 1],
688 z);
689
690 for (iz = 0; iz < model->nz; iz++)
691 zd[iz] = model->z[ilon][ilat + 1][iz];
692 iz = locate_irr(zd, model->nz, z);
693 double p01 = LIN(model->z[ilon][ilat + 1][iz], model->p[ilon][ilat + 1][iz],
694 model->z[ilon][ilat + 1][iz + 1],
695 model->p[ilon][ilat + 1][iz + 1], z);
696 double t01 = LIN(model->z[ilon][ilat + 1][iz], model->t[ilon][ilat + 1][iz],
697 model->z[ilon][ilat + 1][iz + 1],
698 model->t[ilon][ilat + 1][iz + 1],
699 z);
700
701 for (iz = 0; iz < model->nz; iz++)
702 zd[iz] = model->z[ilon + 1][ilat][iz];
703 iz = locate_irr(zd, model->nz, z);
704 double p10 = LIN(model->z[ilon + 1][ilat][iz], model->p[ilon + 1][ilat][iz],
705 model->z[ilon + 1][ilat][iz + 1],
706 model->p[ilon + 1][ilat][iz + 1], z);
707 double t10 = LIN(model->z[ilon + 1][ilat][iz], model->t[ilon + 1][ilat][iz],
708 model->z[ilon + 1][ilat][iz + 1],
709 model->t[ilon + 1][ilat][iz + 1],
710 z);
711
712 for (iz = 0; iz < model->nz; iz++)
713 zd[iz] = model->z[ilon + 1][ilat + 1][iz];
714 iz = locate_irr(zd, model->nz, z);
715 double p11 =
716 LIN(model->z[ilon + 1][ilat + 1][iz], model->p[ilon + 1][ilat + 1][iz],
717 model->z[ilon + 1][ilat + 1][iz + 1],
718 model->p[ilon + 1][ilat + 1][iz + 1], z);
719 double t11 =
720 LIN(model->z[ilon + 1][ilat + 1][iz], model->t[ilon + 1][ilat + 1][iz],
721 model->z[ilon + 1][ilat + 1][iz + 1],
722 model->t[ilon + 1][ilat + 1][iz + 1], z);
723
724 /* Interpolate horizontally... */
725 p00 = LIN(model->lon[ilon], p00, model->lon[ilon + 1], p10, lon);
726 p11 = LIN(model->lon[ilon], p01, model->lon[ilon + 1], p11, lon);
727 *p = LIN(model->lat[ilat], p00, model->lat[ilat + 1], p11, lat);
728
729 t00 = LIN(model->lon[ilon], t00, model->lon[ilon + 1], t10, lon);
730 t11 = LIN(model->lon[ilon], t01, model->lon[ilon + 1], t11, lon);
731 *t = LIN(model->lat[ilat], t00, model->lat[ilat + 1], t11, lat);
732}
733
734/************************************************************************/
735
736void smooth(
737 model_t *model) {
738
739 static float hp[NLON][NLAT], ht[NLON][NLAT], hz[NLON][NLAT];
740
741 static double wx[10], wy[10];
742
743 const int dlon = 3, dlat = 3;
744
745 /* Set weights... */
746 const double dy = RE * M_PI / 180. * fabs(model->lat[1] - model->lat[0]);
747 for (int ilat = 0; ilat <= dlat; ilat++)
748 wy[ilat] = exp(-0.5 * POW2(ilat * dy * 2.35482 / 20.));
749
750 /* Loop over height levels... */
751 for (int iz = 0; iz < model->nz; iz++) {
752
753 /* Write info... */
754 printf("Smoothing level %d / %d ...\n", iz + 1, model->nz);
755
756 /* Copy data... */
757 for (int ilon = 0; ilon < model->nlon; ilon++)
758 for (int ilat = 0; ilat < model->nlat; ilat++) {
759 hp[ilon][ilat] = model->p[ilon][ilat][iz];
760 ht[ilon][ilat] = model->t[ilon][ilat][iz];
761 hz[ilon][ilat] = model->z[ilon][ilat][iz];
762 }
763
764 /* Loop over latitudes... */
765 for (int ilat = 0; ilat < model->nlat; ilat++) {
766
767 /* Set weights... */
768 const double dx =
769 RE * M_PI / 180. * cos(model->lat[ilat] * M_PI / 180.) *
770 fabs(model->lon[1] - model->lon[0]);
771 for (int ilon = 0; ilon <= dlon; ilon++)
772 wx[ilon] = exp(-0.5 * POW2(ilon * dx * 2.35482 / 20.));
773
774 /* Loop over longitudes... */
775 for (int ilon = 0; ilon < model->nlon; ilon++) {
776 float wsum = 0;
777 model->p[ilon][ilat][iz] = 0;
778 model->t[ilon][ilat][iz] = 0;
779 model->z[ilon][ilat][iz] = 0;
780 for (int ilon2 = GSL_MAX(ilon - dlon, 0);
781 ilon2 <= GSL_MIN(ilon + dlon, model->nlon - 1); ilon2++)
782 for (int ilat2 = GSL_MAX(ilat - dlat, 0);
783 ilat2 <= GSL_MIN(ilat + dlat, model->nlat - 1); ilat2++) {
784 float w = (float) (wx[abs(ilon2 - ilon)] * wy[abs(ilat2 - ilat)]);
785 model->p[ilon][ilat][iz] += w * hp[ilon2][ilat2];
786 model->t[ilon][ilat][iz] += w * ht[ilon2][ilat2];
787 model->z[ilon][ilat][iz] += w * hz[ilon2][ilat2];
788 wsum += w;
789 }
790 model->p[ilon][ilat][iz] /= wsum;
791 model->t[ilon][ilat][iz] /= wsum;
792 model->z[ilon][ilat][iz] /= wsum;
793 }
794 }
795 }
796}
797
798/************************************************************************/
799
800void write_nc(
801 char *filename,
802 wave_t *wave) {
803
804 static double help[WX * WY];
805
806 int ncid, dimid[10], lon_id, lat_id, bt_id, pt_id, var_id;
807
808 /* Create netCDF file... */
809 NC(nc_create(filename, NC_CLOBBER, &ncid));
810
811 /* Set dimensions... */
812 NC(nc_def_dim(ncid, "NTRACK", (size_t) wave->ny, &dimid[0]));
813 NC(nc_def_dim(ncid, "NXTRACK", (size_t) wave->nx, &dimid[1]));
814
815 /* Add variables... */
816 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 2, dimid, &lon_id));
817 add_att(ncid, lon_id, "deg", "footprint longitude");
818 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 2, dimid, &lat_id));
819 add_att(ncid, lat_id, "deg", "footprint latitude");
820 NC(nc_def_var(ncid, "bt", NC_FLOAT, 2, dimid, &bt_id));
821 add_att(ncid, bt_id, "K", "brightness temperature");
822 NC(nc_def_var(ncid, "bt_pt", NC_FLOAT, 2, dimid, &pt_id));
823 add_att(ncid, pt_id, "K", "brightness temperature perturbation");
824 NC(nc_def_var(ncid, "bt_var", NC_FLOAT, 2, dimid, &var_id));
825 add_att(ncid, var_id, "K^2", "brightness temperature variance");
826
827 /* Leave define mode... */
828 NC(nc_enddef(ncid));
829
830 /* Write data... */
831 for (int ix = 0; ix < wave->nx; ix++)
832 for (int iy = 0; iy < wave->ny; iy++)
833 help[iy * wave->nx + ix] = wave->lon[ix][iy];
834 NC(nc_put_var_double(ncid, lon_id, help));
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->lat[ix][iy];
838 NC(nc_put_var_double(ncid, lat_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->temp[ix][iy];
842 NC(nc_put_var_double(ncid, bt_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->pt[ix][iy];
846 NC(nc_put_var_double(ncid, pt_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->var[ix][iy];
850 NC(nc_put_var_double(ncid, var_id, help));
851
852 /* Close file... */
853 NC(nc_close(ncid));
854}
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:615
write_nc
void write_nc(char *filename, wave_t *wave)
Write wave struct to netCDF file.
Definition: issifm.c:800
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:736
NLAT
#define NLAT
Maximum number of model latitudes.
Definition: issifm.c:39
NLON
#define NLON
Maximum number of model longitudes.
Definition: issifm.c:36
read_ctl
void read_ctl(int argc, char *argv[], ctl_t *ctl)
Read forward model control parameters.
Definition: jurassic.c:4547
locate_irr
int locate_irr(const double *xx, const int n, const double x)
Find array index for irregular grid.
Definition: jurassic.c:4103
formod
void formod(const ctl_t *ctl, atm_t *atm, obs_t *obs)
Determine ray paths and compute radiative transfer.
Definition: jurassic.c:3026
cart2geo
void cart2geo(const double *x, double *z, double *lon, double *lat)
Convert Cartesian coordinates to geolocation.
Definition: jurassic.c:108
scan_ctl
double scan_ctl(int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
Definition: jurassic.c:5114
read_shape
void read_shape(const char *filename, double *x, double *y, int *n)
Read shape function.
Definition: jurassic.c:4907
climatology
void climatology(const ctl_t *ctl, atm_t *atm)
Interpolate climatological data.
Definition: jurassic.c:123
geo2cart
void geo2cart(const double z, const double lon, const double lat, double *x)
Convert geolocation to Cartesian coordinates.
Definition: jurassic.c:3500
kernel
void kernel(ctl_t *ctl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute Jacobians.
Definition: jurassic.c:4004
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:383
RE
#define RE
Mean radius of Earth [km].
Definition: jurassic.h:299
POW2
#define POW2(x)
Compute x^2.
Definition: jurassic.h:187
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
NSHAPE
#define NSHAPE
Maximum number of shape function grid points.
Definition: jurassic.h:408
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory.
Definition: jurassic.h:121
NP
#define NP
Maximum number of atmospheric data points.
Definition: jurassic.h:363
LOG
#define LOG(level,...)
Print log message.
Definition: jurassic.h:221
LIN
#define LIN(x0, y0, x1, y1, x)
Compute linear interpolation.
Definition: jurassic.h:164
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
atm_t
Atmospheric data.
Definition: jurassic.h:488
atm_t::time
double time[NP]
Time (seconds since 2000-01-01T00:00Z).
Definition: jurassic.h:494
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: jurassic.h:503
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: jurassic.h:500
atm_t::t
double t[NP]
Temperature [K].
Definition: jurassic.h:509
atm_t::np
int np
Number of data points.
Definition: jurassic.h:491
atm_t::z
double z[NP]
Altitude [km].
Definition: jurassic.h:497
atm_t::p
double p[NP]
Pressure [hPa].
Definition: jurassic.h:506
ctl_t
Forward model control parameters.
Definition: jurassic.h:541
ctl_t::nd
int nd
Number of radiance channels.
Definition: jurassic.h:550
ctl_t::write_bbt
int write_bbt
Use brightness temperature instead of radiance (0=no, 1=yes).
Definition: jurassic.h:658
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
obs_t
Observation geometry and radiance data.
Definition: jurassic.h:734
obs_t::rad
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
Definition: jurassic.h:773
obs_t::vpz
double vpz[NR]
View point altitude [km].
Definition: jurassic.h:752
obs_t::vplat
double vplat[NR]
View point latitude [deg].
Definition: jurassic.h:758
obs_t::obslon
double obslon[NR]
Observer longitude [deg].
Definition: jurassic.h:746
obs_t::obslat
double obslat[NR]
Observer latitude [deg].
Definition: jurassic.h:749
obs_t::obsz
double obsz[NR]
Observer altitude [km].
Definition: jurassic.h:743
obs_t::vplon
double vplon[NR]
View point longitude [deg].
Definition: jurassic.h:755
obs_t::nr
int nr
Number of ray paths.
Definition: jurassic.h:737
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 Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4