GPS Code Collection
libgps.c
Go to the documentation of this file.
1/*
2 This file is part of the GPS Code Collection.
3
4 the GPS 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 GPS 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 GPS Code Collection. If not, see
16 <http://www.gnu.org/licenses/>.
17
18 Copyright (C) 2019-2025 Forschungszentrum Juelich GmbH
19*/
20
26#include "libgps.h"
27
28/*****************************************************************************/
29
30void add_var(
31 int ncid,
32 const char *varname,
33 const char *unit,
34 const char *longname,
35 int type,
36 int dimid[],
37 int *varid,
38 int ndims) {
39
40 double dp = GSL_NAN;
41
42 /* Define variable... */
43 NC(nc_def_var(ncid, varname, type, ndims, dimid, varid));
44
45 /* Set long name... */
46 NC(nc_put_att_text(ncid, *varid, "long_name", strlen(longname), longname));
47
48 /* Set units... */
49 NC(nc_put_att_text(ncid, *varid, "units", strlen(unit), unit));
50
51 /* Set fill value... */
52 NC(nc_put_att_double(ncid, *varid, "_FillValue", type, 1, &dp));
53}
54
55/*****************************************************************************/
56
57void detrend_met(
58 gps_t *gps,
59 char *metbase,
60 double dt_met) {
61
62 met_t *met0, *met1;
63
64 double t;
65
66 /* Allocate... */
67 ALLOC(met0, met_t, 1);
68 ALLOC(met1, met_t, 1);
69
70 /* Loop over profiles... */
71 for (int ids = 0; ids < gps->nds; ids++) {
72
73 /* Loop over altitudes... */
74 for (int iz = 0; iz < gps->nz[ids]; iz++) {
75
76 /* Get meteorological data... */
77 get_met(metbase, dt_met, gps->time[ids], met0, met1);
78
79 /* Interpolate meteorological data... */
80 intpol_met_time(met0, met1, gps->time[ids], gps->p[ids][iz],
81 gps->lon[ids][iz], gps->lat[ids][iz], &t);
82
83 /* Set perturbation... */
84 gps->pt[ids][iz] = gps->t[ids][iz] - t;
85 }
86 }
87
88 /* Free... */
89 free(met0);
90 free(met1);
91}
92
93/*****************************************************************************/
94
95void gauss(
96 gps_t *gps,
97 double dx,
98 double dy) {
99
100 /* Loop over profiles... */
101 for (int ids = 0; ids < gps->nds; ids++) {
102
103 /* Initialize... */
104 double wsum = 0;
105 for (int iz = 0; iz < gps->nz[ids]; iz++)
106 gps->pt[ids][iz] = 0;
107
108 /* Calculate lon-lat standard deviations... */
109 const double dlat = dx * 180. / (M_PI * RE) / 2.3548;
110 const double dlon = dy * 180. / 2.3548
111 / (M_PI * RE * cos(gps->lat[ids][gps->nz[ids] / 2] * M_PI / 180.));
112
113 /* Calculate mean temperature... */
114 for (int ids2 = 0; ids2 < gps->nds; ids2++) {
115 const double w = exp(-0.5 * gsl_pow_2((gps->lon[ids][gps->nz[ids] / 2]
116 - gps->lon[ids2][gps->nz[ids2] / 2]) / dlon)
117 - 0.5 * gsl_pow_2((gps->lat[ids][gps->nz[ids] / 2]
118 - gps->lat[ids2][gps->nz[ids2] / 2]) / dlat));
119 wsum += w;
120 for (int iz = 0; iz < gps->nz[ids]; iz++)
121 gps->pt[ids][iz] += w * gps->t[ids2][iz];
122 }
123
124 /* Normalize... */
125 if (wsum > 0)
126 for (int iz = 0; iz < gps->nz[ids]; iz++)
127 gps->pt[ids][iz] = gps->t[ids][iz] - gps->pt[ids][iz] / wsum;
128 }
129}
130
131/*****************************************************************************/
132
133void grid_gps(
134 gps_t *gps,
135 double zmin,
136 double zmax,
137 int nz) {
138
139 double lat[NZ], lon[NZ], p[NZ], pt[NZ], t[NZ], wv[NZ], z[NZ];
140
141 /* Check number of altitudes... */
142 if (nz > NZ)
143 ERRMSG("Too many altitudes!");
144
145 /* Loop over profiles... */
146 for (int ids = 0; ids < gps->nds; ids++) {
147
148 /* Loop over altitudes... */
149 for (int iz = 0; iz < nz; iz++) {
150
151 /* Set altitude... */
152 z[iz] = LIN(0.0, zmin, nz - 1.0, zmax, (double) iz);
153
154 /* Get index... */
155 const int iz2 = locate_irr(gps->z[ids], gps->nz[ids], z[iz]);
156
157 /* Interpolate... */
158 lon[iz] = LIN(gps->z[ids][iz2], gps->lon[ids][iz2],
159 gps->z[ids][iz2 + 1], gps->lon[ids][iz2 + 1], z[iz]);
160 lat[iz] = LIN(gps->z[ids][iz2], gps->lat[ids][iz2],
161 gps->z[ids][iz2 + 1], gps->lat[ids][iz2 + 1], z[iz]);
162 p[iz] = LIN(gps->z[ids][iz2], gps->p[ids][iz2],
163 gps->z[ids][iz2 + 1], gps->p[ids][iz2 + 1], z[iz]);
164 t[iz] = LIN(gps->z[ids][iz2], gps->t[ids][iz2],
165 gps->z[ids][iz2 + 1], gps->t[ids][iz2 + 1], z[iz]);
166 wv[iz] = LIN(gps->z[ids][iz2], gps->wv[ids][iz2],
167 gps->z[ids][iz2 + 1], gps->wv[ids][iz2 + 1], z[iz]);
168 pt[iz] = LIN(gps->z[ids][iz2], gps->pt[ids][iz2],
169 gps->z[ids][iz2 + 1], gps->pt[ids][iz2 + 1], z[iz]);
170 }
171
172 /* Copy data... */
173 gps->nz[ids] = nz;
174 for (int iz = 0; iz < nz; iz++) {
175 gps->z[ids][iz] = z[iz];
176 gps->lon[ids][iz] = lon[iz];
177 gps->lat[ids][iz] = lat[iz];
178 gps->p[ids][iz] = p[iz];
179 gps->t[ids][iz] = t[iz];
180 gps->wv[ids][iz] = wv[iz];
181 gps->pt[ids][iz] = pt[iz];
182 }
183 }
184}
185
186/*****************************************************************************/
187
188void get_met(
189 char *metbase,
190 double dt_met,
191 double t,
192 met_t *met0,
193 met_t *met1) {
194
195 char filename[LEN];
196
197 static int init;
198
199 /* Init... */
200 if (!init) {
201 init = 1;
202
203 get_met_help(t, -1, metbase, dt_met, filename);
204 read_met(filename, met0);
205
206 get_met_help(t + 1.0, 1, metbase, dt_met, filename);
207 read_met(filename, met1);
208 }
209
210 /* Read new data... */
211 if (t > met1->time) {
212 memcpy(met0, met1, sizeof(met_t));
213 get_met_help(t, 1, metbase, dt_met, filename);
214 read_met(filename, met1);
215 }
216}
217
218/*****************************************************************************/
219
220void get_met_help(
221 double t,
222 int direct,
223 char *metbase,
224 double dt_met,
225 char *filename) {
226
227 double t6, r;
228
229 int year, mon, day, hour, min, sec;
230
231 /* Round time to fixed intervals... */
232 if (direct == -1)
233 t6 = floor(t / dt_met) * dt_met;
234 else
235 t6 = ceil(t / dt_met) * dt_met;
236
237 /* Decode time... */
238 jsec2time(t6, &year, &mon, &day, &hour, &min, &sec, &r);
239
240 /* Set filename... */
241 sprintf(filename, "%s_%d_%02d_%02d_%02d.nc", metbase, year, mon, day, hour);
242}
243
244/*****************************************************************************/
245
246void intpol_met_3d(
247 float array[EX][EY][EP],
248 int ip,
249 int ix,
250 int iy,
251 double wp,
252 double wx,
253 double wy,
254 double *var) {
255
256 /* Interpolate vertically... */
257 double aux00 = wp * (array[ix][iy][ip] - array[ix][iy][ip + 1])
258 + array[ix][iy][ip + 1];
259 double aux01 = wp * (array[ix][iy + 1][ip] - array[ix][iy + 1][ip + 1])
260 + array[ix][iy + 1][ip + 1];
261 double aux10 = wp * (array[ix + 1][iy][ip] - array[ix + 1][iy][ip + 1])
262 + array[ix + 1][iy][ip + 1];
263 double aux11 = wp * (array[ix + 1][iy + 1][ip] - array[ix + 1][iy + 1][ip + 1])
264 + array[ix + 1][iy + 1][ip + 1];
265
266 /* Interpolate horizontally... */
267 aux00 = wy * (aux00 - aux01) + aux01;
268 aux11 = wy * (aux10 - aux11) + aux11;
269 *var = wx * (aux00 - aux11) + aux11;
270}
271
272/*****************************************************************************/
273
274void intpol_met_space(
275 met_t *met,
276 double p,
277 double lon,
278 double lat,
279 double *t) {
280
281 /* Check longitude... */
282 if (met->lon[met->nx - 1] > 180 && lon < 0)
283 lon += 360;
284
285 /* Get indices... */
286 const int ip = locate_irr(met->p, met->np, p);
287 const int ix = locate_reg(met->lon, met->nx, lon);
288 const int iy = locate_reg(met->lat, met->ny, lat);
289
290 /* Get weights... */
291 const double wp = (met->p[ip + 1] - p) / (met->p[ip + 1] - met->p[ip]);
292 const double wx = (met->lon[ix + 1] - lon) / (met->lon[ix + 1] - met->lon[ix]);
293 const double wy = (met->lat[iy + 1] - lat) / (met->lat[iy + 1] - met->lat[iy]);
294
295 /* Interpolate... */
296 intpol_met_3d(met->t, ip, ix, iy, wp, wx, wy, t);
297}
298
299/*****************************************************************************/
300
301void intpol_met_time(
302 met_t *met0,
303 met_t *met1,
304 double ts,
305 double p,
306 double lon,
307 double lat,
308 double *t) {
309
310 double t0, t1;
311
312 /* Spatial interpolation... */
313 intpol_met_space(met0, p, lon, lat, &t0);
314 intpol_met_space(met1, p, lon, lat, &t1);
315
316 /* Get weighting factor... */
317 const double wt = (met1->time - ts) / (met1->time - met0->time);
318
319 /* Interpolate... */
320 *t = wt * (t0 - t1) + t1;
321}
322
323/*****************************************************************************/
324
325void hamming_low_pass(
326 gps_t *gps,
327 double dz) {
328
329 double ham[NZ];
330
331 /* Loop over profiles... */
332 for (int ids = 0; ids < gps->nds; ids++) {
333
334 /* Calculate Hamming window coefficients... */
335 int nham = (int) (dz / fabs((gps->z[ids][0] - gps->z[ids][gps->nz[ids] - 1])
336 / (gps->nz[ids] - 1.0)) + 0.5);
337 nham = GSL_MAX(GSL_MIN(nham, NZ), 2);
338 for (int iham = 0; iham < nham; iham++)
339 ham[iham] = 0.54 + 0.46 * cos(M_PI * iham / (nham - 1.0));
340
341 /* Loop over altitudes... */
342 for (int iz = 0; iz < gps->nz[ids]; iz++) {
343
344 /* Initialize... */
345 gps->pt[ids][iz] = ham[0] * gps->t[ids][iz];
346 double wsum = ham[0];
347
348 /* Loop over filter window... */
349 for (int iham = 1; iham < nham; iham++) {
350
351 /* Check array range... */
352 if (iz - iham < 0 || iz + iham >= gps->nz[ids])
353 continue;
354
355 /* Check temperature value... */
356 if (!gsl_finite(gps->t[ids][iz - iham]) ||
357 !gsl_finite(gps->t[ids][iz + iham]))
358 continue;
359
360 /* Check for tropopause... */
361 if (gsl_finite(gps->th[ids]) && gps->th[ids] > 0)
362 if ((gps->z[ids][iz] >= gps->th[ids]
363 && gps->z[ids][iz - iham] < gps->th[ids])
364 || (gps->z[ids][iz] <= gps->th[ids]
365 && gps->z[ids][iz + iham] > gps->th[ids]))
366 continue;
367
368 /* Apply Hamming filter... */
369 gps->pt[ids][iz]
370 += ham[iham] * (gps->t[ids][iz - iham] + gps->t[ids][iz + iham]);
371 wsum += 2 * ham[iham];
372 }
373
374 /* Calculate perturbation... */
375 gps->pt[ids][iz] = gps->t[ids][iz] - gps->pt[ids][iz] / wsum;
376 }
377 }
378}
379
380/*****************************************************************************/
381
382void hamming_high_pass(
383 gps_t *gps,
384 double dz) {
385
386 double ham[NZ], pt[NZ];
387
388 /* Loop over profiles... */
389 for (int ids = 0; ids < gps->nds; ids++) {
390
391 /* Calculate Hamming window coefficients... */
392 int nham = (int) (dz / fabs((gps->z[ids][0] - gps->z[ids][gps->nz[ids] - 1])
393 / (gps->nz[ids] - 1.0)) + 0.5);
394 nham = GSL_MAX(GSL_MIN(nham, NZ), 2);
395 for (int iham = 0; iham < nham; iham++)
396 ham[iham] = 0.54 + 0.46 * cos(M_PI * iham / (nham - 1.0));
397
398 /* Loop over altitudes... */
399 for (int iz = 0; iz < gps->nz[ids]; iz++) {
400
401 /* Initialize... */
402 pt[iz] = ham[0] * gps->pt[ids][iz];
403 double wsum = ham[0];
404
405 /* Loop over filter window... */
406 for (int iham = 1; iham < nham; iham++) {
407
408 /* Check array range... */
409 if (iz - iham < 0 || iz + iham >= gps->nz[ids])
410 continue;
411
412 /* Check temperature value... */
413 if (!gsl_finite(gps->t[ids][iz - iham]) ||
414 !gsl_finite(gps->t[ids][iz + iham]))
415 continue;
416
417 /* Apply Hamming filter... */
418 pt[iz]
419 += ham[iham] * (gps->pt[ids][iz - iham] + gps->pt[ids][iz + iham]);
420 wsum += 2 * ham[iham];
421 }
422
423 /* Normalize... */
424 pt[iz] /= wsum;
425 }
426
427 /* Set perturbation... */
428 for (int iz = 0; iz < gps->nz[ids]; iz++)
429 gps->pt[ids][iz] = pt[iz];
430 }
431}
432
433/*****************************************************************************/
434
435void poly(
436 gps_t *gps,
437 int dim,
438 double zmin,
439 double zmax) {
440
441 double bg[NZ];
442
443 /* Loop over profiles... */
444 for (int ids = 0; ids < gps->nds; ids++) {
445
446 /* Set profile... */
447 for (int iz = 0; iz < gps->nz[ids]; iz++)
448 bg[iz] = gps->pt[ids][iz];
449
450 /* Polynomial interpolation... */
451 poly_help(gps->z[ids], bg, gps->nz[ids], dim, zmin, zmax);
452
453 /* Remove background... */
454 for (int iz = 0; iz < gps->nz[ids]; iz++)
455 gps->pt[ids][iz] -= bg[iz];
456 }
457}
458
459/*****************************************************************************/
460
461void poly_help(
462 double *xx,
463 double *yy,
464 int n,
465 int dim,
466 double xmin,
467 double xmax) {
468
469 gsl_multifit_linear_workspace *work;
470 gsl_matrix *cov, *X;
471 gsl_vector *c, *x, *y;
472
473 double chisq, xx2[NZ], yy2[NZ];
474
475 size_t i, i2, n2 = 0;
476
477 /* Check for nan... */
478 for (i = 0; i < (size_t) n; i++)
479 if (xx[i] >= xmin && xx[i] <= xmax && gsl_finite(yy[i])) {
480 xx2[n2] = xx[i];
481 yy2[n2] = yy[i];
482 n2++;
483 }
484 if ((int) n2 < dim) {
485 for (i = 0; i < (size_t) n; i++)
486 yy[i] = GSL_NAN;
487 return;
488 }
489
490 /* Allocate... */
491 work = gsl_multifit_linear_alloc((size_t) n2, (size_t) dim);
492 cov = gsl_matrix_alloc((size_t) dim, (size_t) dim);
493 X = gsl_matrix_alloc((size_t) n2, (size_t) dim);
494 c = gsl_vector_alloc((size_t) dim);
495 x = gsl_vector_alloc((size_t) n2);
496 y = gsl_vector_alloc((size_t) n2);
497
498 /* Compute polynomial fit... */
499 for (i = 0; i < (size_t) n2; i++) {
500 gsl_vector_set(x, i, xx2[i]);
501 gsl_vector_set(y, i, yy2[i]);
502 for (i2 = 0; i2 < (size_t) dim; i2++)
503 gsl_matrix_set(X, i, i2, pow(gsl_vector_get(x, i), (double) i2));
504 }
505 gsl_multifit_linear(X, y, c, cov, &chisq, work);
506 for (i = 0; i < (size_t) n; i++)
507 yy[i] = gsl_poly_eval(c->data, (int) dim, xx[i]);
508
509 /* Free... */
510 gsl_multifit_linear_free(work);
511 gsl_matrix_free(cov);
512 gsl_matrix_free(X);
513 gsl_vector_free(c);
514 gsl_vector_free(x);
515 gsl_vector_free(y);
516}
517
518/*****************************************************************************/
519
520void read_gps_prof(
521 char *filename,
522 gps_t *gps) {
523
524 char bad[10];
525
526 double t0, t1, zmin = 1e100, zmax = -1e100;
527
528 int ncid, dimid, varid;
529
530 size_t iz, nz;
531
532 /* Open netCDF file... */
533 printf("Read GPS-RO profile: %s\n", filename);
534 NC(nc_open(filename, NC_NOWRITE, &ncid));
535
536 /* Get dimensions... */
537 NC(nc_inq_dimid(ncid, "MSL_alt", &dimid));
538 NC(nc_inq_dimlen(ncid, dimid, &nz));
539 gps->nz[gps->nds] = (int) nz;
540 if (nz > NZ)
541 ERRMSG("Too many altitudes!");
542
543 /* Check data quality flag... */
544 NC(nc_get_att_text(ncid, NC_GLOBAL, "bad", bad));
545 if (bad[0] != '0') {
546 NC(nc_close(ncid));
547 return;
548 }
549
550 /* Get time... */
551 NC(nc_get_att_double(ncid, NC_GLOBAL, "start_time", &t0));
552 NC(nc_get_att_double(ncid, NC_GLOBAL, "stop_time", &t1));
553 gps->time[gps->nds] = 0.5 * (t0 + t1) - 630720000.0;
554
555 /* Get data... */
556 NC(nc_inq_varid(ncid, "MSL_alt", &varid));
557 NC(nc_get_var_double(ncid, varid, gps->z[gps->nds]));
558 NC(nc_inq_varid(ncid, "Lon", &varid));
559 NC(nc_get_var_double(ncid, varid, gps->lon[gps->nds]));
560 NC(nc_inq_varid(ncid, "Lat", &varid));
561 NC(nc_get_var_double(ncid, varid, gps->lat[gps->nds]));
562 NC(nc_inq_varid(ncid, "Pres", &varid));
563 NC(nc_get_var_double(ncid, varid, gps->p[gps->nds]));
564 NC(nc_inq_varid(ncid, "Temp", &varid));
565 NC(nc_get_var_double(ncid, varid, gps->t[gps->nds]));
566 if (nc_inq_varid(ncid, "Vp", &varid) == NC_NOERR)
567 NC(nc_get_var_double(ncid, varid, gps->wv[gps->nds]));
568
569 /* Check altitude range... */
570 for (iz = 0; iz < nz; iz++)
571 if (gps->p[gps->nds][iz] != -999 && gps->t[gps->nds][iz] != -999) {
572 zmin = GSL_MIN(zmin, gps->z[gps->nds][iz]);
573 zmax = GSL_MAX(zmax, gps->z[gps->nds][iz]);
574 }
575 if (zmin > 5 || zmax < 35) {
576 NC(nc_close(ncid));
577 return;
578 }
579
580 /* Check data... */
581 for (iz = 0; iz < nz; iz++)
582 if (gps->lon[gps->nds][iz] == -999 ||
583 gps->lat[gps->nds][iz] == -999 ||
584 gps->p[gps->nds][iz] == -999 ||
585 gps->t[gps->nds][iz] == -999 || gps->wv[gps->nds][iz] == -999) {
586 gps->lon[gps->nds][iz] = GSL_NAN;
587 gps->lat[gps->nds][iz] = GSL_NAN;
588 gps->p[gps->nds][iz] = GSL_NAN;
589 gps->t[gps->nds][iz] = GSL_NAN;
590 gps->wv[gps->nds][iz] = GSL_NAN;
591 }
592
593 /* Convert temperature... */
594 for (iz = 0; iz < nz; iz++)
595 gps->t[gps->nds][iz] += 273.15;
596
597 /* Convert water vapor... */
598 for (iz = 0; iz < nz; iz++)
599 gps->wv[gps->nds][iz] /= gps->p[gps->nds][iz];
600
601 /* Close file... */
602 NC(nc_close(ncid));
603
604 /* Count profiles... */
605 if ((++gps->nds) >= NDS)
606 ERRMSG("Too many profiles!");
607}
608
609/*****************************************************************************/
610
611void read_gps(
612 char *filename,
613 gps_t *gps) {
614
615 int ids, ncid, dimid, varid;
616
617 size_t start[2], count[2], nds, nz;
618
619 /* Read netCDF file... */
620 printf("Read GPS-RO file: %s\n", filename);
621 NC(nc_open(filename, NC_NOWRITE, &ncid));
622
623 /* Get dimensions... */
624 NC(nc_inq_dimid(ncid, "NDS", &dimid));
625 NC(nc_inq_dimlen(ncid, dimid, &nds));
626 gps->nds = (int) nds;
627 if (nds > NDS)
628 ERRMSG("Too many profiles!");
629
630 NC(nc_inq_dimid(ncid, "NZ", &dimid));
631 NC(nc_inq_dimlen(ncid, dimid, &nz));
632 if (nz > NZ)
633 ERRMSG("Too many profiles!");
634
635 /* Loop over profiles... */
636 for (ids = 0; ids < gps->nds; ids++) {
637
638 /* Set profile index... */
639 start[0] = (size_t) ids;
640 count[0] = 1;
641 start[1] = 0;
642 count[1] = nz;
643
644 /* Set number of altitudes... */
645 gps->nz[ids] = (int) nz;
646
647 /* Read data... */
648 NC(nc_inq_varid(ncid, "time", &varid));
649 NC(nc_get_vara_double(ncid, varid, start, count, &gps->time[ids]));
650
651 NC(nc_inq_varid(ncid, "z", &varid));
652 NC(nc_get_vara_double(ncid, varid, start, count, gps->z[ids]));
653
654 NC(nc_inq_varid(ncid, "lon", &varid));
655 NC(nc_get_vara_double(ncid, varid, start, count, gps->lon[ids]));
656
657 NC(nc_inq_varid(ncid, "lat", &varid));
658 NC(nc_get_vara_double(ncid, varid, start, count, gps->lat[ids]));
659
660 NC(nc_inq_varid(ncid, "p", &varid));
661 NC(nc_get_vara_double(ncid, varid, start, count, gps->p[ids]));
662
663 NC(nc_inq_varid(ncid, "t", &varid));
664 NC(nc_get_vara_double(ncid, varid, start, count, gps->t[ids]));
665
666 NC(nc_inq_varid(ncid, "wv", &varid));
667 NC(nc_get_vara_double(ncid, varid, start, count, gps->wv[ids]));
668
669 NC(nc_inq_varid(ncid, "pt", &varid));
670 NC(nc_get_vara_double(ncid, varid, start, count, gps->pt[ids]));
671
672 NC(nc_inq_varid(ncid, "th", &varid));
673 NC(nc_get_vara_double(ncid, varid, start, count, &gps->th[ids]));
674 }
675
676 /* Close file... */
677 NC(nc_close(ncid));
678}
679
680/*****************************************************************************/
681
682void read_met(
683 char *filename,
684 met_t *met) {
685
686 char tstr[10];
687
688 int ip, dimid, ncid, varid, year, mon, day, hour;
689
690 size_t np, nx, ny;
691
692 /* Write info... */
693 printf("Read meteorological data: %s\n", filename);
694
695 /* Get time from filename... */
696 sprintf(tstr, "%.4s", &filename[strlen(filename) - 16]);
697 year = atoi(tstr);
698 sprintf(tstr, "%.2s", &filename[strlen(filename) - 11]);
699 mon = atoi(tstr);
700 sprintf(tstr, "%.2s", &filename[strlen(filename) - 8]);
701 day = atoi(tstr);
702 sprintf(tstr, "%.2s", &filename[strlen(filename) - 5]);
703 hour = atoi(tstr);
704 time2jsec(year, mon, day, hour, 0, 0, 0, &met->time);
705
706 /* Open netCDF file... */
707 NC(nc_open(filename, NC_NOWRITE, &ncid));
708
709 /* Get dimensions... */
710 NC(nc_inq_dimid(ncid, "lon", &dimid));
711 NC(nc_inq_dimlen(ncid, dimid, &nx));
712 if (nx > EX)
713 ERRMSG("Too many longitudes!");
714
715 NC(nc_inq_dimid(ncid, "lat", &dimid));
716 NC(nc_inq_dimlen(ncid, dimid, &ny));
717 if (ny > EY)
718 ERRMSG("Too many latitudes!");
719
720 NC(nc_inq_dimid(ncid, "lev", &dimid));
721 NC(nc_inq_dimlen(ncid, dimid, &np));
722 if (np > EP)
723 ERRMSG("Too many levels!");
724
725 /* Store dimensions... */
726 met->np = (int) np;
727 met->nx = (int) nx;
728 met->ny = (int) ny;
729
730 /* Get horizontal grid... */
731 NC(nc_inq_varid(ncid, "lon", &varid));
732 NC(nc_get_var_double(ncid, varid, met->lon));
733 NC(nc_inq_varid(ncid, "lat", &varid));
734 NC(nc_get_var_double(ncid, varid, met->lat));
735
736 /* Read meteorological data... */
737 read_met_help(ncid, "t", "T", met, met->t, 1.0);
738
739 /* Read pressure levels from file... */
740 NC(nc_inq_varid(ncid, "lev", &varid));
741 NC(nc_get_var_double(ncid, varid, met->p));
742 for (ip = 0; ip < met->np; ip++)
743 met->p[ip] /= 100.;
744
745 /* Extrapolate data for lower boundary... */
746 read_met_extrapolate(met);
747
748 /* Check ordering of pressure levels... */
749 for (ip = 1; ip < met->np; ip++)
750 if (met->p[ip - 1] < met->p[ip])
751 ERRMSG("Pressure levels must be descending!");
752
753 /* Create periodic boundary conditions... */
754 read_met_periodic(met);
755
756 /* Close file... */
757 NC(nc_close(ncid));
758}
759
760/*****************************************************************************/
761
762void read_met_extrapolate(
763 met_t *met) {
764
765 int ip0;
766
767 /* Loop over columns... */
768 for (int ix = 0; ix < met->nx; ix++)
769 for (int iy = 0; iy < met->ny; iy++) {
770
771 /* Find lowest valid data point... */
772 for (ip0 = met->np - 1; ip0 >= 0; ip0--)
773 if (!gsl_finite(met->t[ix][iy][ip0]))
774 break;
775
776 /* Extrapolate... */
777 for (int ip = ip0; ip >= 0; ip--)
778 met->t[ix][iy][ip] = met->t[ix][iy][ip + 1];
779 }
780}
781
782/*****************************************************************************/
783
784void read_met_help(
785 int ncid,
786 char *varname,
787 char *varname2,
788 met_t *met,
789 float dest[EX][EY][EP],
790 float scl) {
791
792 static float help[EX * EY * EP];
793
794 int n = 0, varid;
795
796 /* Check if variable exists... */
797 if (nc_inq_varid(ncid, varname, &varid) != NC_NOERR)
798 if (nc_inq_varid(ncid, varname2, &varid) != NC_NOERR)
799 return;
800
801 /* Read data... */
802 NC(nc_get_var_float(ncid, varid, help));
803
804 /* Copy and check data... */
805 for (int ip = 0; ip < met->np; ip++)
806 for (int iy = 0; iy < met->ny; iy++)
807 for (int ix = 0; ix < met->nx; ix++) {
808 dest[ix][iy][ip] = scl * help[n++];
809 if (fabs(dest[ix][iy][ip] / scl) > 1e14)
810 dest[ix][iy][ip] = GSL_NAN;
811 }
812}
813
814/*****************************************************************************/
815
816void read_met_periodic(
817 met_t *met) {
818
819 /* Check longitudes... */
820 if (!(fabs(met->lon[met->nx - 1] - met->lon[0]
821 + met->lon[1] - met->lon[0] - 360) < 0.01))
822 return;
823
824 /* Increase longitude counter... */
825 if ((++met->nx) > EX)
826 ERRMSG("Cannot create periodic boundary conditions!");
827
828 /* Set longitude... */
829 met->lon[met->nx - 1] = met->lon[met->nx - 2] + met->lon[1] - met->lon[0];
830
831 /* Loop over latitudes and pressure levels... */
832 for (int iy = 0; iy < met->ny; iy++)
833 for (int ip = 0; ip < met->np; ip++)
834 met->t[met->nx - 1][iy][ip] = met->t[0][iy][ip];
835}
836
837/*****************************************************************************/
838
839void spline(
840 const double *x,
841 const double *y,
842 const int n,
843 const double *x2,
844 double *y2,
845 const int n2,
846 const int method) {
847
848 /* Cubic spline interpolation... */
849 if (method == 1) {
850
851 /* Allocate... */
852 gsl_interp_accel *acc = gsl_interp_accel_alloc();
853 gsl_spline *s = gsl_spline_alloc(gsl_interp_cspline, (size_t) n);
854
855 /* Interpolate profile... */
856 gsl_spline_init(s, x, y, (size_t) n);
857 for (int i = 0; i < n2; i++)
858 if (x2[i] <= x[0])
859 y2[i] = y[0];
860 else if (x2[i] >= x[n - 1])
861 y2[i] = y[n - 1];
862 else
863 y2[i] = gsl_spline_eval(s, x2[i], acc);
864
865 /* Free... */
866 gsl_spline_free(s);
867 gsl_interp_accel_free(acc);
868 }
869
870 /* Linear interpolation... */
871 else {
872 for (int i = 0; i < n2; i++)
873 if (x2[i] <= x[0])
874 y2[i] = y[0];
875 else if (x2[i] >= x[n - 1])
876 y2[i] = y[n - 1];
877 else {
878 int idx = locate_irr(x, n, x2[i]);
879 y2[i] = LIN(x[idx], y[idx], x[idx + 1], y[idx + 1], x2[i]);
880 }
881 }
882}
883
884/*****************************************************************************/
885
886void tropopause(
887 gps_t *gps) {
888
889 int iz, iz2, okay;
890
891 /* Loop over profiles... */
892 for (int ids = 0; ids < gps->nds; ids++) {
893
894 /* Set default value... */
895 gps->th[ids] = GSL_NAN;
896
897 /* Set minimum altitude... */
898 const double zmin =
899 8 - 4 * fabs(cos((90 - gps->lat[ids][gps->nz[ids] / 2]) * M_PI / 180));
900
901 /* Search tropopause (WMO definition)... */
902 for (iz = 0; iz < gps->nz[ids]; iz++)
903 if (gps->z[ids][iz] >= zmin && gps->z[ids][iz] <= 20.0) {
904 okay = 1;
905 for (iz2 = iz + 1; iz2 < gps->nz[ids]; iz2++)
906 if (gps->z[ids][iz2] - gps->z[ids][iz] <= 2.0)
907 if (!gsl_finite(gps->t[ids][iz]) ||
908 !gsl_finite(gps->t[ids][iz2]) ||
909 (gps->t[ids][iz2] - gps->t[ids][iz])
910 / (gps->z[ids][iz2] - gps->z[ids][iz]) < -2.0)
911 okay = 0;
912 if (okay) {
913 gps->th[ids] = gps->z[ids][iz];
914 break;
915 }
916 }
917 }
918}
919
920/*****************************************************************************/
921
922void tropopause_spline(
923 gps_t *gps,
924 int met_tropo) {
925
926 /* Loop over data sets... */
927#pragma omp parallel for default(shared)
928 for (int ids = 0; ids < gps->nds; ids++) {
929
930 /* Init... */
931 gps->tp[ids] = NAN;
932 gps->th[ids] = NAN;
933 gps->tt[ids] = NAN;
934 gps->tq[ids] = NAN;
935 gps->tlon[ids] = NAN;
936 gps->tlat[ids] = NAN;
937
938 /* Get vertical profiles... */
939 int nz = 0;
940 double h[NZ], t[NZ], z[NZ], q[NZ], lon[NZ], lat[NZ];
941 for (int iz = 0; iz < gps->nz[ids]; iz++)
942 if (gsl_finite(gps->p[ids][iz]) && gsl_finite(gps->t[ids][iz])
943 && gsl_finite(gps->z[ids][iz])
944 && gps->z[ids][iz] >= 4.0 && gps->z[ids][iz] <= 24.0) {
945 h[nz] = gps->z[ids][iz];
946 t[nz] = gps->t[ids][iz];
947 z[nz] = Z(gps->p[ids][iz]);
948 q[nz] = gps->wv[ids][iz];
949 lon[nz] = gps->lon[ids][iz];
950 lat[nz] = gps->lat[ids][iz];
951 if (nz > 0 && z[nz] <= z[nz - 1])
952 ERRMSG("Profiles must be ascending!");
953 if ((++nz) >= NZ)
954 ERRMSG("Too many height levels!");
955 }
956 if (z[0] > 4.5 || z[nz - 1] < 23.5)
957 WARN("Vertical profile is incomplete!");
958
959 /* Set grid for spline interpolation... */
960 double h2[200], p2[200], t2[200], z2[200], q2[200];
961 for (int iz = 0; iz <= 190; iz++) {
962 z2[iz] = 4.5 + 0.1 * iz;
963 p2[iz] = P(z2[iz]);
964 }
965
966 /* Interpolate temperature and geopotential height profiles... */
967 spline(z, t, nz, z2, t2, 191, 1);
968 spline(z, h, nz, z2, h2, 191, 1);
969 spline(z, q, nz, z2, q2, 191, 1);
970
971 /* Use cold point... */
972 if (met_tropo == 2) {
973
974 /* Find minimum... */
975 int iz = (int) gsl_stats_min_index(t2, 1, 171);
976 if (iz > 0 && iz < 170) {
977 gps->tp[ids] = p2[iz];
978 gps->th[ids] = h2[iz];
979 gps->tt[ids] = t2[iz];
980 gps->tq[ids] = q2[iz];
981 }
982 }
983
984 /* Use WMO definition... */
985 else if (met_tropo == 3 || met_tropo == 4) {
986
987 /* Find 1st tropopause... */
988 int iz;
989 for (iz = 0; iz <= 170; iz++) {
990 int found = 1;
991 for (int iz2 = iz + 1; iz2 <= iz + 20; iz2++)
992 if (LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]) > 2.0) {
993 found = 0;
994 break;
995 }
996 if (found) {
997 if (iz > 0 && iz < 170) {
998 gps->tp[ids] = p2[iz];
999 gps->th[ids] = h2[iz];
1000 gps->tt[ids] = t2[iz];
1001 gps->tq[ids] = q2[iz];
1002 }
1003 break;
1004 }
1005 }
1006
1007 /* Find 2nd tropopause... */
1008 if (met_tropo == 4) {
1009
1010 /* Init... */
1011 gps->tp[ids] = NAN;
1012 gps->th[ids] = NAN;
1013 gps->tt[ids] = NAN;
1014 gps->tq[ids] = NAN;
1015
1016 /* Check layers... */
1017 for (; iz <= 170; iz++) {
1018 int found = 1;
1019 for (int iz2 = iz + 1; iz2 <= iz + 10; iz2++)
1020 if (LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]) < 3.0) {
1021 found = 0;
1022 break;
1023 }
1024 if (found)
1025 break;
1026 }
1027 for (; iz <= 170; iz++) {
1028 int found = 1;
1029 for (int iz2 = iz + 1; iz2 <= iz + 20; iz2++)
1030 if (LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]) > 2.0) {
1031 found = 0;
1032 break;
1033 }
1034 if (found) {
1035 if (iz > 0 && iz < 170) {
1036 gps->tp[ids] = p2[iz];
1037 gps->th[ids] = h2[iz];
1038 gps->tt[ids] = t2[iz];
1039 gps->tq[ids] = q2[iz];
1040 }
1041 break;
1042 }
1043 }
1044 }
1045 }
1046
1047 /* Find tropopause longitude and latitude... */
1048 if (gsl_finite(gps->th[ids]))
1049 for (int iz = 0; iz < nz - 1; iz++)
1050 if (gps->th[ids] >= h[iz] && gps->th[ids] < h[iz + 1]) {
1051 gps->tlon[ids] = lon[iz];
1052 gps->tlat[ids] = lat[iz];
1053 break;
1054 }
1055 }
1056}
1057
1058/*****************************************************************************/
1059
1060void write_gps(
1061 char *filename,
1062 gps_t *gps) {
1063
1064 static double help[NDS * NZ];
1065
1066 int ids, iz, ncid, dimid[2], time_id, z_id, lon_id, lat_id, p_id, t_id,
1067 pt_id, wv_id, th_id, nzmax = 0;
1068
1069 /* Create netCDF file... */
1070 printf("Write GPS-RO file: %s\n", filename);
1071 NC(nc_create(filename, NC_CLOBBER, &ncid));
1072
1073 /* Set dimensions... */
1074 NC(nc_def_dim(ncid, "NDS", (size_t) gps->nds, &dimid[0]));
1075 for (ids = 0; ids < gps->nds; ids++)
1076 nzmax = GSL_MAX(nzmax, gps->nz[ids]);
1077 NC(nc_def_dim(ncid, "NZ", (size_t) nzmax, &dimid[1]));
1078
1079 /* Add variables... */
1080 add_var(ncid, "time", "s", "time (seconds since 2000-01-01T00:00Z)",
1081 NC_DOUBLE, dimid, &time_id, 1);
1082 add_var(ncid, "z", "km", "altitude", NC_FLOAT, dimid, &z_id, 2);
1083 add_var(ncid, "lon", "deg", "longitude", NC_FLOAT, dimid, &lon_id, 2);
1084 add_var(ncid, "lat", "deg", "latitude", NC_FLOAT, dimid, &lat_id, 2);
1085 add_var(ncid, "p", "hPa", "pressure", NC_FLOAT, dimid, &p_id, 2);
1086 add_var(ncid, "t", "K", "temperature", NC_FLOAT, dimid, &t_id, 2);
1087 add_var(ncid, "wv", "ppv", "water vapor volume mixing ratio",
1088 NC_FLOAT, dimid, &wv_id, 2);
1089 add_var(ncid, "pt", "K", "temperature perturbation",
1090 NC_FLOAT, dimid, &pt_id, 2);
1091 add_var(ncid, "th", "km", "tropopause height", NC_FLOAT, dimid, &th_id, 1);
1092
1093 /* Leave define mode... */
1094 NC(nc_enddef(ncid));
1095
1096 /* Write data... */
1097 NC(nc_put_var_double(ncid, time_id, gps->time));
1098 NC(nc_put_var_double(ncid, th_id, gps->th));
1099 for (ids = 0; ids < gps->nds; ids++)
1100 for (iz = 0; iz < gps->nz[ids]; iz++)
1101 help[ids * nzmax + iz] = gps->z[ids][iz];
1102 NC(nc_put_var_double(ncid, z_id, help));
1103 for (ids = 0; ids < gps->nds; ids++)
1104 for (iz = 0; iz < gps->nz[ids]; iz++)
1105 help[ids * nzmax + iz] = gps->lon[ids][iz];
1106 NC(nc_put_var_double(ncid, lon_id, help));
1107 for (ids = 0; ids < gps->nds; ids++)
1108 for (iz = 0; iz < gps->nz[ids]; iz++)
1109 help[ids * nzmax + iz] = gps->lat[ids][iz];
1110 NC(nc_put_var_double(ncid, lat_id, help));
1111 for (ids = 0; ids < gps->nds; ids++)
1112 for (iz = 0; iz < gps->nz[ids]; iz++)
1113 help[ids * nzmax + iz] = gps->p[ids][iz];
1114 NC(nc_put_var_double(ncid, p_id, help));
1115 for (ids = 0; ids < gps->nds; ids++)
1116 for (iz = 0; iz < gps->nz[ids]; iz++)
1117 help[ids * nzmax + iz] = gps->t[ids][iz];
1118 NC(nc_put_var_double(ncid, t_id, help));
1119 for (ids = 0; ids < gps->nds; ids++)
1120 for (iz = 0; iz < gps->nz[ids]; iz++)
1121 help[ids * nzmax + iz] = gps->wv[ids][iz];
1122 NC(nc_put_var_double(ncid, wv_id, help));
1123 for (ids = 0; ids < gps->nds; ids++)
1124 for (iz = 0; iz < gps->nz[ids]; iz++)
1125 help[ids * nzmax + iz] = gps->pt[ids][iz];
1126 NC(nc_put_var_double(ncid, pt_id, help));
1127
1128 /* Close file... */
1129 NC(nc_close(ncid));
1130}
read_met_extrapolate
void read_met_extrapolate(met_t *met)
Extrapolate meteorological data at lower boundary.
Definition: libgps.c:762
intpol_met_space
void intpol_met_space(met_t *met, double p, double lon, double lat, double *t)
Spatial interpolation of meteorological data.
Definition: libgps.c:274
write_gps
void write_gps(char *filename, gps_t *gps)
Write GPS-RO data file.
Definition: libgps.c:1060
read_met_periodic
void read_met_periodic(met_t *met)
Create meteorological data with periodic boundary conditions.
Definition: libgps.c:816
gauss
void gauss(gps_t *gps, double dx, double dy)
Calculate horizontal Gaussian mean to extract perturbations.
Definition: libgps.c:95
read_met
void read_met(char *filename, met_t *met)
Read meteorological data file.
Definition: libgps.c:682
intpol_met_time
void intpol_met_time(met_t *met0, met_t *met1, double ts, double p, double lon, double lat, double *t)
Temporal interpolation of meteorological data.
Definition: libgps.c:301
spline
void spline(const double *x, const double *y, const int n, const double *x2, double *y2, const int n2, const int method)
Performs spline interpolation or linear interpolation.
Definition: libgps.c:839
poly
void poly(gps_t *gps, int dim, double zmin, double zmax)
Remove polynomial fit from perturbation profile.
Definition: libgps.c:435
grid_gps
void grid_gps(gps_t *gps, double zmin, double zmax, int nz)
Interpolate GPS data to regular altitude grid.
Definition: libgps.c:133
intpol_met_3d
void intpol_met_3d(float array[EX][EY][EP], int ip, int ix, int iy, double wp, double wx, double wy, double *var)
Linear interpolation of 3-D meteorological data.
Definition: libgps.c:246
add_var
void add_var(int ncid, const char *varname, const char *unit, const char *longname, int type, int dimid[], int *varid, int ndims)
Add variable to netCDF file.
Definition: libgps.c:30
get_met
void get_met(char *metbase, double dt_met, double t, met_t *met0, met_t *met1)
Get meteorological data for given timestep.
Definition: libgps.c:188
read_gps
void read_gps(char *filename, gps_t *gps)
Read GPS-RO data file.
Definition: libgps.c:611
hamming_low_pass
void hamming_low_pass(gps_t *gps, double dz)
Apply vertical Hamming filter to extract perturbations.
Definition: libgps.c:325
detrend_met
void detrend_met(gps_t *gps, char *metbase, double dt_met)
Detrending by means of meteo data.
Definition: libgps.c:57
read_gps_prof
void read_gps_prof(char *filename, gps_t *gps)
Read GPS-RO profile.
Definition: libgps.c:520
get_met_help
void get_met_help(double t, int direct, char *metbase, double dt_met, char *filename)
Get meteorological data for timestep.
Definition: libgps.c:220
hamming_high_pass
void hamming_high_pass(gps_t *gps, double dz)
Apply vertical Hamming filter to reduce noise.
Definition: libgps.c:382
tropopause_spline
void tropopause_spline(gps_t *gps, int met_tropo)
Find tropopause height using cubic spline interpolation.
Definition: libgps.c:922
poly_help
void poly_help(double *xx, double *yy, int n, int dim, double xmin, double xmax)
Auxiliary function for polynomial interpolation.
Definition: libgps.c:461
read_met_help
void read_met_help(int ncid, char *varname, char *varname2, met_t *met, float dest[EX][EY][EP], float scl)
Read and convert variable from meteorological data file.
Definition: libgps.c:784
tropopause
void tropopause(gps_t *gps)
Find tropopause height.
Definition: libgps.c:886
libgps.h
GPS Code Collection library declarations.
LAPSE
#define LAPSE(p1, t1, p2, t2)
Calculate lapse rate.
Definition: libgps.h:129
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libgps.h:134
EY
#define EY
Maximum number of latitudes for meteorological data.
Definition: libgps.h:97
Z
#define Z(p)
Convert pressure to altitude.
Definition: libgps.h:145
P
#define P(z)
Compute pressure at given altitude.
Definition: libgps.h:141
NZ
#define NZ
Maximum number of altitudes per GPS-RO profile.
Definition: libgps.h:103
EX
#define EX
Maximum number of longitudes for meteorological data.
Definition: libgps.h:94
NDS
#define NDS
Maximum number of GPS-RO profiles.
Definition: libgps.h:100
EP
#define EP
Maximum number of pressure levels for meteorological data.
Definition: libgps.h:91
gps_t
GPS-RO profile data.
Definition: libgps.h:153
gps_t::time
double time[NDS]
Time (seconds since 2000-01-01T00:00Z).
Definition: libgps.h:162
gps_t::pt
double pt[NDS][NZ]
Temperature perturbation [K].
Definition: libgps.h:183
gps_t::tlon
double tlon[NDS]
Tropopause longitude [deg].
Definition: libgps.h:198
gps_t::tt
double tt[NDS]
Tropopause temperature [K].
Definition: libgps.h:192
gps_t::t
double t[NDS][NZ]
Temperature [K].
Definition: libgps.h:177
gps_t::tp
double tp[NDS]
Tropopause pressure [hPa].
Definition: libgps.h:189
gps_t::wv
double wv[NDS][NZ]
Water vapor volume mixing ratio [ppm].
Definition: libgps.h:180
gps_t::nz
int nz[NDS]
Number of altitudes per profile.
Definition: libgps.h:159
gps_t::lon
double lon[NDS][NZ]
Longitude [deg].
Definition: libgps.h:168
gps_t::nds
int nds
Number of profiles.
Definition: libgps.h:156
gps_t::p
double p[NDS][NZ]
Pressure [hPa].
Definition: libgps.h:174
gps_t::th
double th[NDS]
Tropopause height [km].
Definition: libgps.h:186
gps_t::z
double z[NDS][NZ]
Altitude [km].
Definition: libgps.h:165
gps_t::tlat
double tlat[NDS]
Tropopause latitude [deg].
Definition: libgps.h:201
gps_t::lat
double lat[NDS][NZ]
Latitude [deg].
Definition: libgps.h:171
gps_t::tq
double tq[NDS]
Tropopause water vapor [ppmv].
Definition: libgps.h:195
met_t
Meteorological data.
Definition: libgps.h:206
met_t::nx
int nx
Number of longitudes.
Definition: libgps.h:212
met_t::ny
int ny
Number of latitudes.
Definition: libgps.h:215
met_t::np
int np
Number of pressure levels.
Definition: libgps.h:218
met_t::t
float t[EX][EY][EP]
Temperature [K].
Definition: libgps.h:230
met_t::lon
double lon[EX]
Longitude [deg].
Definition: libgps.h:221
met_t::time
double time
Time [s].
Definition: libgps.h:209
met_t::lat
double lat[EY]
Latitude [deg].
Definition: libgps.h:224
met_t::p
double p[EP]
Pressure [hPa].
Definition: libgps.h:227
Generated on Sat Jul 19 2025 19:02:12 for GPS Code Collection by doxygen 1.9.4