GPS Code Collection
Functions
libgps.c File Reference

GPS Code Collection library definitions. More...

#include "libgps.h"

Go to the source code of this file.

Functions

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. More...
 
void detrend_met (gps_t *gps, char *metbase, double dt_met)
 Detrending by means of meteo data. More...
 
void gauss (gps_t *gps, double dx, double dy)
 Calculate horizontal Gaussian mean to extract perturbations. More...
 
void grid_gps (gps_t *gps, double zmin, double zmax, int nz)
 Interpolate GPS data to regular altitude grid. More...
 
void get_met (char *metbase, double dt_met, double t, met_t *met0, met_t *met1)
 Get meteorological data for given timestep. More...
 
void get_met_help (double t, int direct, char *metbase, double dt_met, char *filename)
 Get meteorological data for timestep. More...
 
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. More...
 
void intpol_met_space (met_t *met, double p, double lon, double lat, double *t)
 Spatial interpolation of meteorological data. More...
 
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. More...
 
void hamming_low_pass (gps_t *gps, double dz)
 Apply vertical Hamming filter to extract perturbations. More...
 
void hamming_high_pass (gps_t *gps, double dz)
 Apply vertical Hamming filter to reduce noise. More...
 
void poly (gps_t *gps, int dim, double zmin, double zmax)
 Remove polynomial fit from perturbation profile. More...
 
void poly_help (double *xx, double *yy, int n, int dim, double xmin, double xmax)
 Auxiliary function for polynomial interpolation. More...
 
void read_gps_prof (char *filename, gps_t *gps, double prof_zmin_max, double prof_zmax_min)
 Read GPS-RO profile. More...
 
void read_gps (char *filename, gps_t *gps)
 Read GPS-RO data file. More...
 
void read_met (char *filename, met_t *met)
 Read meteorological data file. More...
 
void read_met_extrapolate (met_t *met)
 Extrapolate meteorological data at lower boundary. More...
 
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. More...
 
void read_met_periodic (met_t *met)
 Create meteorological data with periodic boundary conditions. More...
 
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. More...
 
void tropopause (gps_t *gps)
 Find tropopause height. More...
 
void tropopause_spline (gps_t *gps, int met_tropo)
 Find tropopause height using cubic spline interpolation. More...
 
void write_gps (char *filename, gps_t *gps)
 Write GPS-RO data file. More...
 

Detailed Description

GPS Code Collection library definitions.

Definition in file libgps.c.

Function Documentation

◆ 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 at line 30 of file libgps.c.

38 {
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}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libgps.h:134

◆ detrend_met()

void detrend_met ( gps_t *  gps,
char *  metbase,
double  dt_met 
)

Detrending by means of meteo data.

Definition at line 57 of file libgps.c.

60 {
61
62 met_t *met0, *met1;
63
64 /* Allocate... */
65 ALLOC(met0, met_t, 1);
66 ALLOC(met1, met_t, 1);
67
68 /* Loop over profiles... */
69 for (int ids = 0; ids < gps->nds; ids++) {
70
71 /* Loop over altitudes... */
72 for (int iz = 0; iz < gps->nz[ids]; iz++) {
73
74 /* Get meteorological data... */
75 get_met(metbase, dt_met, gps->time[ids], met0, met1);
76
77 /* Interpolate meteorological data... */
78 double t;
79 intpol_met_time(met0, met1, gps->time[ids], gps->p[ids][iz],
80 gps->lon[ids][iz], gps->lat[ids][iz], &t);
81
82 /* Set perturbation... */
83 gps->pt[ids][iz] = gps->t[ids][iz] - t;
84 }
85 }
86
87 /* Free... */
88 free(met0);
89 free(met1);
90}
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:307
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:191
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::t
double t[NDS][NZ]
Temperature [K].
Definition: libgps.h:177
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::lat
double lat[NDS][NZ]
Latitude [deg].
Definition: libgps.h:171
met_t
Meteorological data.
Definition: libgps.h:206
Here is the call graph for this function:

◆ gauss()

void gauss ( gps_t *  gps,
double  dx,
double  dy 
)

Calculate horizontal Gaussian mean to extract perturbations.

Definition at line 94 of file libgps.c.

97 {
98
99 /* Loop over profiles... */
100 for (int ids = 0; ids < gps->nds; ids++) {
101
102 /* Initialize... */
103 double wsum = 0;
104 for (int iz = 0; iz < gps->nz[ids]; iz++)
105 gps->pt[ids][iz] = 0;
106
107 /* Calculate lon-lat standard deviations... */
108 const double dlat = dx * 180. / (M_PI * RE) / 2.3548;
109 const double dlon = dy * 180. / 2.3548
110 / (M_PI * RE * cos(gps->lat[ids][gps->nz[ids] / 2] * M_PI / 180.));
111
112 /* Calculate mean temperature... */
113 for (int ids2 = 0; ids2 < gps->nds; ids2++) {
114 const double w = exp(-0.5 * gsl_pow_2((gps->lon[ids][gps->nz[ids] / 2]
115 -
116 gps->lon[ids2][gps->nz[ids2] /
117 2]) / dlon)
118 - 0.5 * gsl_pow_2((gps->lat[ids][gps->nz[ids] / 2]
119 -
120 gps->lat[ids2][gps->nz[ids2] /
121 2]) / dlat));
122 wsum += w;
123 for (int iz = 0; iz < gps->nz[ids]; iz++)
124 gps->pt[ids][iz] += w * gps->t[ids2][iz];
125 }
126
127 /* Normalize... */
128 if (wsum > 0)
129 for (int iz = 0; iz < gps->nz[ids]; iz++)
130 gps->pt[ids][iz] = gps->t[ids][iz] - gps->pt[ids][iz] / wsum;
131 }
132}

◆ grid_gps()

void grid_gps ( gps_t *  gps,
double  zmin,
double  zmax,
int  nz 
)

Interpolate GPS data to regular altitude grid.

Definition at line 136 of file libgps.c.

140 {
141
142 double lat[NZ], lon[NZ], p[NZ], pt[NZ], t[NZ], wv[NZ], z[NZ];
143
144 /* Check number of altitudes... */
145 if (nz > NZ)
146 ERRMSG("Too many altitudes!");
147
148 /* Loop over profiles... */
149 for (int ids = 0; ids < gps->nds; ids++) {
150
151 /* Loop over altitudes... */
152 for (int iz = 0; iz < nz; iz++) {
153
154 /* Set altitude... */
155 z[iz] = LIN(0.0, zmin, nz - 1.0, zmax, (double) iz);
156
157 /* Get index... */
158 const int iz2 = locate_irr(gps->z[ids], gps->nz[ids], z[iz]);
159
160 /* Interpolate... */
161 lon[iz] = LIN(gps->z[ids][iz2], gps->lon[ids][iz2],
162 gps->z[ids][iz2 + 1], gps->lon[ids][iz2 + 1], z[iz]);
163 lat[iz] = LIN(gps->z[ids][iz2], gps->lat[ids][iz2],
164 gps->z[ids][iz2 + 1], gps->lat[ids][iz2 + 1], z[iz]);
165 p[iz] = LIN(gps->z[ids][iz2], gps->p[ids][iz2],
166 gps->z[ids][iz2 + 1], gps->p[ids][iz2 + 1], z[iz]);
167 t[iz] = LIN(gps->z[ids][iz2], gps->t[ids][iz2],
168 gps->z[ids][iz2 + 1], gps->t[ids][iz2 + 1], z[iz]);
169 wv[iz] = LIN(gps->z[ids][iz2], gps->wv[ids][iz2],
170 gps->z[ids][iz2 + 1], gps->wv[ids][iz2 + 1], z[iz]);
171 pt[iz] = LIN(gps->z[ids][iz2], gps->pt[ids][iz2],
172 gps->z[ids][iz2 + 1], gps->pt[ids][iz2 + 1], z[iz]);
173 }
174
175 /* Copy data... */
176 gps->nz[ids] = nz;
177 for (int iz = 0; iz < nz; iz++) {
178 gps->z[ids][iz] = z[iz];
179 gps->lon[ids][iz] = lon[iz];
180 gps->lat[ids][iz] = lat[iz];
181 gps->p[ids][iz] = p[iz];
182 gps->t[ids][iz] = t[iz];
183 gps->wv[ids][iz] = wv[iz];
184 gps->pt[ids][iz] = pt[iz];
185 }
186 }
187}
NZ
#define NZ
Maximum number of altitudes per GPS-RO profile.
Definition: libgps.h:103
gps_t::wv
double wv[NDS][NZ]
Water vapor volume mixing ratio [ppm].
Definition: libgps.h:180
gps_t::z
double z[NDS][NZ]
Altitude [km].
Definition: libgps.h:165

◆ 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 at line 191 of file libgps.c.

196 {
197
198 char filename[LEN];
199
200 static int init;
201
202 /* Init... */
203 if (!init) {
204 init = 1;
205
206 get_met_help(t, -1, metbase, dt_met, filename);
207 read_met(filename, met0);
208
209 get_met_help(t + 1.0, 1, metbase, dt_met, filename);
210 read_met(filename, met1);
211 }
212
213 /* Read new data... */
214 if (t > met1->time) {
215 memcpy(met0, met1, sizeof(met_t));
216 get_met_help(t, 1, metbase, dt_met, filename);
217 read_met(filename, met1);
218 }
219}
read_met
void read_met(char *filename, met_t *met)
Read meteorological data file.
Definition: libgps.c:705
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:223
met_t::time
double time
Time [s].
Definition: libgps.h:209
Here is the call graph for this function:

◆ get_met_help()

void get_met_help ( double  t,
int  direct,
char *  metbase,
double  dt_met,
char *  filename 
)

Get meteorological data for timestep.

Definition at line 223 of file libgps.c.

228 {
229
230 double t6, r;
231
232 int year, mon, day, hour, min, sec;
233
234 /* Round time to fixed intervals... */
235 if (direct == -1)
236 t6 = floor(t / dt_met) * dt_met;
237 else
238 t6 = ceil(t / dt_met) * dt_met;
239
240 /* Decode time... */
241 jsec2time(t6, &year, &mon, &day, &hour, &min, &sec, &r);
242
243 /* Set filename... */
244 sprintf(filename, "%s_%d_%02d_%02d_%02d.nc", metbase, year, mon, day, hour);
245}

◆ 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 at line 249 of file libgps.c.

257 {
258
259 /* Interpolate vertically... */
260 double aux00 = wp * (array[ix][iy][ip] - array[ix][iy][ip + 1])
261 + array[ix][iy][ip + 1];
262 double aux01 = wp * (array[ix][iy + 1][ip] - array[ix][iy + 1][ip + 1])
263 + array[ix][iy + 1][ip + 1];
264 double aux10 = wp * (array[ix + 1][iy][ip] - array[ix + 1][iy][ip + 1])
265 + array[ix + 1][iy][ip + 1];
266 double aux11 =
267 wp * (array[ix + 1][iy + 1][ip] - array[ix + 1][iy + 1][ip + 1])
268 + array[ix + 1][iy + 1][ip + 1];
269
270 /* Interpolate horizontally... */
271 aux00 = wy * (aux00 - aux01) + aux01;
272 aux11 = wy * (aux10 - aux11) + aux11;
273 *var = wx * (aux00 - aux11) + aux11;
274}

◆ intpol_met_space()

void intpol_met_space ( met_t *  met,
double  p,
double  lon,
double  lat,
double *  t 
)

Spatial interpolation of meteorological data.

Definition at line 278 of file libgps.c.

283 {
284
285 /* Check longitude... */
286 if (met->lon[met->nx - 1] > 180 && lon < 0)
287 lon += 360;
288
289 /* Get indices... */
290 const int ip = locate_irr(met->p, met->np, p);
291 const int ix = locate_reg(met->lon, met->nx, lon);
292 const int iy = locate_reg(met->lat, met->ny, lat);
293
294 /* Get weights... */
295 const double wp = (met->p[ip + 1] - p) / (met->p[ip + 1] - met->p[ip]);
296 const double wx =
297 (met->lon[ix + 1] - lon) / (met->lon[ix + 1] - met->lon[ix]);
298 const double wy =
299 (met->lat[iy + 1] - lat) / (met->lat[iy + 1] - met->lat[iy]);
300
301 /* Interpolate... */
302 intpol_met_3d(met->t, ip, ix, iy, wp, wx, wy, t);
303}
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:249
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::lat
double lat[EY]
Latitude [deg].
Definition: libgps.h:224
met_t::p
double p[EP]
Pressure [hPa].
Definition: libgps.h:227
Here is the call graph for this function:

◆ 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 at line 307 of file libgps.c.

314 {
315
316 double t0, t1;
317
318 /* Spatial interpolation... */
319 intpol_met_space(met0, p, lon, lat, &t0);
320 intpol_met_space(met1, p, lon, lat, &t1);
321
322 /* Get weighting factor... */
323 const double wt = (met1->time - ts) / (met1->time - met0->time);
324
325 /* Interpolate... */
326 *t = wt * (t0 - t1) + t1;
327}
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:278
Here is the call graph for this function:

◆ hamming_low_pass()

void hamming_low_pass ( gps_t *  gps,
double  dz 
)

Apply vertical Hamming filter to extract perturbations.

Definition at line 331 of file libgps.c.

333 {
334
335 double ham[NZ];
336
337 /* Loop over profiles... */
338 for (int ids = 0; ids < gps->nds; ids++) {
339
340 /* Calculate Hamming window coefficients... */
341 int nham =
342 (int) (dz / fabs((gps->z[ids][0] - gps->z[ids][gps->nz[ids] - 1])
343 / (gps->nz[ids] - 1.0)) + 0.5);
344 nham = GSL_MAX(GSL_MIN(nham, NZ), 2);
345 for (int iham = 0; iham < nham; iham++)
346 ham[iham] = 0.54 + 0.46 * cos(M_PI * iham / (nham - 1.0));
347
348 /* Loop over altitudes... */
349 for (int iz = 0; iz < gps->nz[ids]; iz++) {
350
351 /* Initialize... */
352 gps->pt[ids][iz] = ham[0] * gps->t[ids][iz];
353 double wsum = ham[0];
354
355 /* Loop over filter window... */
356 for (int iham = 1; iham < nham; iham++) {
357
358 /* Check array range... */
359 if (iz - iham < 0 || iz + iham >= gps->nz[ids])
360 continue;
361
362 /* Check temperature value... */
363 if (!gsl_finite(gps->t[ids][iz - iham]) ||
364 !gsl_finite(gps->t[ids][iz + iham]))
365 continue;
366
367 /* Check for tropopause... */
368 if (gsl_finite(gps->th[ids]) && gps->th[ids] > 0)
369 if ((gps->z[ids][iz] >= gps->th[ids]
370 && gps->z[ids][iz - iham] < gps->th[ids])
371 || (gps->z[ids][iz] <= gps->th[ids]
372 && gps->z[ids][iz + iham] > gps->th[ids]))
373 continue;
374
375 /* Apply Hamming filter... */
376 gps->pt[ids][iz]
377 += ham[iham] * (gps->t[ids][iz - iham] + gps->t[ids][iz + iham]);
378 wsum += 2 * ham[iham];
379 }
380
381 /* Calculate perturbation... */
382 gps->pt[ids][iz] = gps->t[ids][iz] - gps->pt[ids][iz] / wsum;
383 }
384 }
385}
gps_t::th
double th[NDS]
Tropopause height [km].
Definition: libgps.h:186

◆ hamming_high_pass()

void hamming_high_pass ( gps_t *  gps,
double  dz 
)

Apply vertical Hamming filter to reduce noise.

Definition at line 389 of file libgps.c.

391 {
392
393 double ham[NZ], pt[NZ];
394
395 /* Loop over profiles... */
396 for (int ids = 0; ids < gps->nds; ids++) {
397
398 /* Calculate Hamming window coefficients... */
399 int nham =
400 (int) (dz / fabs((gps->z[ids][0] - gps->z[ids][gps->nz[ids] - 1])
401 / (gps->nz[ids] - 1.0)) + 0.5);
402 nham = GSL_MAX(GSL_MIN(nham, NZ), 2);
403 for (int iham = 0; iham < nham; iham++)
404 ham[iham] = 0.54 + 0.46 * cos(M_PI * iham / (nham - 1.0));
405
406 /* Loop over altitudes... */
407 for (int iz = 0; iz < gps->nz[ids]; iz++) {
408
409 /* Initialize... */
410 pt[iz] = ham[0] * gps->pt[ids][iz];
411 double wsum = ham[0];
412
413 /* Loop over filter window... */
414 for (int iham = 1; iham < nham; iham++) {
415
416 /* Check array range... */
417 if (iz - iham < 0 || iz + iham >= gps->nz[ids])
418 continue;
419
420 /* Check temperature value... */
421 if (!gsl_finite(gps->t[ids][iz - iham]) ||
422 !gsl_finite(gps->t[ids][iz + iham]))
423 continue;
424
425 /* Apply Hamming filter... */
426 pt[iz]
427 += ham[iham] * (gps->pt[ids][iz - iham] + gps->pt[ids][iz + iham]);
428 wsum += 2 * ham[iham];
429 }
430
431 /* Normalize... */
432 pt[iz] /= wsum;
433 }
434
435 /* Set perturbation... */
436 for (int iz = 0; iz < gps->nz[ids]; iz++)
437 gps->pt[ids][iz] = pt[iz];
438 }
439}

◆ poly()

void poly ( gps_t *  gps,
int  dim,
double  zmin,
double  zmax 
)

Remove polynomial fit from perturbation profile.

Definition at line 443 of file libgps.c.

447 {
448
449 double bg[NZ];
450
451 /* Loop over profiles... */
452 for (int ids = 0; ids < gps->nds; ids++) {
453
454 /* Set profile... */
455 for (int iz = 0; iz < gps->nz[ids]; iz++)
456 bg[iz] = gps->pt[ids][iz];
457
458 /* Polynomial interpolation... */
459 poly_help(gps->z[ids], bg, gps->nz[ids], dim, zmin, zmax);
460
461 /* Remove background... */
462 for (int iz = 0; iz < gps->nz[ids]; iz++)
463 gps->pt[ids][iz] -= bg[iz];
464 }
465}
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:469
Here is the call graph for this function:

◆ poly_help()

void poly_help ( double *  xx,
double *  yy,
int  n,
int  dim,
double  xmin,
double  xmax 
)

Auxiliary function for polynomial interpolation.

Definition at line 469 of file libgps.c.

475 {
476
477 gsl_multifit_linear_workspace *work;
478 gsl_matrix *cov, *X;
479 gsl_vector *c, *x, *y;
480
481 double chisq, xx2[NZ], yy2[NZ];
482
483 size_t n2 = 0;
484
485 /* Check for nan... */
486 for (size_t i = 0; i < (size_t) n; i++)
487 if (xx[i] >= xmin && xx[i] <= xmax && gsl_finite(yy[i])) {
488 xx2[n2] = xx[i];
489 yy2[n2] = yy[i];
490 n2++;
491 }
492 if ((int) n2 < dim) {
493 for (size_t i = 0; i < (size_t) n; i++)
494 yy[i] = GSL_NAN;
495 return;
496 }
497
498 /* Allocate... */
499 work = gsl_multifit_linear_alloc((size_t) n2, (size_t) dim);
500 cov = gsl_matrix_alloc((size_t) dim, (size_t) dim);
501 X = gsl_matrix_alloc((size_t) n2, (size_t) dim);
502 c = gsl_vector_alloc((size_t) dim);
503 x = gsl_vector_alloc((size_t) n2);
504 y = gsl_vector_alloc((size_t) n2);
505
506 /* Compute polynomial fit... */
507 for (size_t i = 0; i < (size_t) n2; i++) {
508 gsl_vector_set(x, i, xx2[i]);
509 gsl_vector_set(y, i, yy2[i]);
510 for (size_t i2 = 0; i2 < (size_t) dim; i2++)
511 gsl_matrix_set(X, i, i2, pow(gsl_vector_get(x, i), (double) i2));
512 }
513 gsl_multifit_linear(X, y, c, cov, &chisq, work);
514 for (size_t i = 0; i < (size_t) n; i++)
515 yy[i] = gsl_poly_eval(c->data, (int) dim, xx[i]);
516
517 /* Free... */
518 gsl_multifit_linear_free(work);
519 gsl_matrix_free(cov);
520 gsl_matrix_free(X);
521 gsl_vector_free(c);
522 gsl_vector_free(x);
523 gsl_vector_free(y);
524}

◆ read_gps_prof()

void read_gps_prof ( char *  filename,
gps_t *  gps,
double  prof_zmin_max,
double  prof_zmax_min 
)

Read GPS-RO profile.

Definition at line 528 of file libgps.c.

532 {
533
534 char bad[10];
535
536 double t0, t1, zmin = 1e100, zmax = -1e100;
537
538 int ncid, dimid, dimid2, varid;
539
540 size_t nz;
541
542 /* Open netCDF file... */
543 printf("Read GPS-RO profile: %s\n", filename);
544 NC(nc_open(filename, NC_NOWRITE, &ncid));
545
546 /* Get dimensions... */
547 NC(nc_inq_dimid(ncid, "MSL_alt", &dimid));
548 NC(nc_inq_dimlen(ncid, dimid, &nz));
549 gps->nz[gps->nds] = (int) nz;
550 if (nz > NZ)
551 ERRMSG("Too many altitudes!");
552
553 /* Check data quality flag... */
554 if (nc_get_att_text(ncid, NC_GLOBAL, "bad", bad) == NC_NOERR)
555 if (bad[0] != '0') {
556 NC(nc_close(ncid));
557 return;
558 }
559
560 /* Get time... */
561 if (nc_get_att_double(ncid, NC_GLOBAL, "start_time", &t0) == NC_NOERR
562 && nc_get_att_double(ncid, NC_GLOBAL, "stop_time", &t1) == NC_NOERR)
563 gps->time[gps->nds] = 0.5 * (t0 + t1) - 630720000.0;
564 else {
565 NC(nc_inq_varid(ncid, "Time", &varid));
566 NC(nc_get_var_double(ncid, varid, &gps->time[gps->nds]));
567 gps->time[gps->nds] -= 630720000.0;
568 }
569
570 /* Get data... */
571 NC(nc_inq_varid(ncid, "MSL_alt", &varid));
572 NC(nc_get_var_double(ncid, varid, gps->z[gps->nds]));
573 NC(nc_inq_varid(ncid, "Lon", &varid));
574 NC(nc_get_var_double(ncid, varid, gps->lon[gps->nds]));
575 NC(nc_inq_var(ncid, varid, NULL, NULL, NULL, &dimid2, NULL));
576 if (dimid2 != dimid)
577 for (size_t iz = 1; iz < nz; iz++)
578 gps->lon[gps->nds][iz] = gps->lon[gps->nds][0];
579 NC(nc_inq_varid(ncid, "Lat", &varid));
580 NC(nc_get_var_double(ncid, varid, gps->lat[gps->nds]));
581 NC(nc_inq_var(ncid, varid, NULL, NULL, NULL, &dimid2, NULL));
582 if (dimid2 != dimid)
583 for (size_t iz = 1; iz < nz; iz++)
584 gps->lat[gps->nds][iz] = gps->lat[gps->nds][0];
585 NC(nc_inq_varid(ncid, "Pres", &varid));
586 NC(nc_get_var_double(ncid, varid, gps->p[gps->nds]));
587 NC(nc_inq_varid(ncid, "Temp", &varid));
588 NC(nc_get_var_double(ncid, varid, gps->t[gps->nds]));
589 if (nc_inq_varid(ncid, "Vp", &varid) == NC_NOERR)
590 NC(nc_get_var_double(ncid, varid, gps->wv[gps->nds]));
591
592 /* Check altitude range... */
593 for (size_t iz = 0; iz < nz; iz++)
594 if (gps->p[gps->nds][iz] != -999 && gps->t[gps->nds][iz] != -999) {
595 zmin = GSL_MIN(zmin, gps->z[gps->nds][iz]);
596 zmax = GSL_MAX(zmax, gps->z[gps->nds][iz]);
597 }
598 if (zmin > prof_zmin_max || zmax < prof_zmax_min) {
599 NC(nc_close(ncid));
600 return;
601 }
602
603 /* Check data... */
604 for (size_t iz = 0; iz < nz; iz++)
605 if (gps->lon[gps->nds][iz] == -999 ||
606 gps->lat[gps->nds][iz] == -999 ||
607 gps->p[gps->nds][iz] == -999 ||
608 gps->t[gps->nds][iz] == -999 || gps->wv[gps->nds][iz] == -999) {
609 gps->lon[gps->nds][iz] = GSL_NAN;
610 gps->lat[gps->nds][iz] = GSL_NAN;
611 gps->p[gps->nds][iz] = GSL_NAN;
612 gps->t[gps->nds][iz] = GSL_NAN;
613 gps->wv[gps->nds][iz] = GSL_NAN;
614 }
615
616 /* Convert temperature... */
617 for (size_t iz = 0; iz < nz; iz++)
618 gps->t[gps->nds][iz] += 273.15;
619
620 /* Convert water vapor... */
621 for (size_t iz = 0; iz < nz; iz++)
622 gps->wv[gps->nds][iz] /= gps->p[gps->nds][iz];
623
624 /* Close file... */
625 NC(nc_close(ncid));
626
627 /* Count profiles... */
628 if ((++gps->nds) >= NDS)
629 ERRMSG("Too many profiles!");
630}
NDS
#define NDS
Maximum number of GPS-RO profiles.
Definition: libgps.h:100

◆ read_gps()

void read_gps ( char *  filename,
gps_t *  gps 
)

Read GPS-RO data file.

Definition at line 634 of file libgps.c.

636 {
637
638 int ids, ncid, dimid, varid;
639
640 size_t start[2], count[2], nds, nz;
641
642 /* Read netCDF file... */
643 printf("Read GPS-RO file: %s\n", filename);
644 NC(nc_open(filename, NC_NOWRITE, &ncid));
645
646 /* Get dimensions... */
647 NC(nc_inq_dimid(ncid, "NDS", &dimid));
648 NC(nc_inq_dimlen(ncid, dimid, &nds));
649 gps->nds = (int) nds;
650 if (nds > NDS)
651 ERRMSG("Too many profiles!");
652
653 NC(nc_inq_dimid(ncid, "NZ", &dimid));
654 NC(nc_inq_dimlen(ncid, dimid, &nz));
655 if (nz > NZ)
656 ERRMSG("Too many profiles!");
657
658 /* Loop over profiles... */
659 for (ids = 0; ids < gps->nds; ids++) {
660
661 /* Set profile index... */
662 start[0] = (size_t) ids;
663 count[0] = 1;
664 start[1] = 0;
665 count[1] = nz;
666
667 /* Set number of altitudes... */
668 gps->nz[ids] = (int) nz;
669
670 /* Read data... */
671 NC(nc_inq_varid(ncid, "time", &varid));
672 NC(nc_get_vara_double(ncid, varid, start, count, &gps->time[ids]));
673
674 NC(nc_inq_varid(ncid, "z", &varid));
675 NC(nc_get_vara_double(ncid, varid, start, count, gps->z[ids]));
676
677 NC(nc_inq_varid(ncid, "lon", &varid));
678 NC(nc_get_vara_double(ncid, varid, start, count, gps->lon[ids]));
679
680 NC(nc_inq_varid(ncid, "lat", &varid));
681 NC(nc_get_vara_double(ncid, varid, start, count, gps->lat[ids]));
682
683 NC(nc_inq_varid(ncid, "p", &varid));
684 NC(nc_get_vara_double(ncid, varid, start, count, gps->p[ids]));
685
686 NC(nc_inq_varid(ncid, "t", &varid));
687 NC(nc_get_vara_double(ncid, varid, start, count, gps->t[ids]));
688
689 NC(nc_inq_varid(ncid, "wv", &varid));
690 NC(nc_get_vara_double(ncid, varid, start, count, gps->wv[ids]));
691
692 NC(nc_inq_varid(ncid, "pt", &varid));
693 NC(nc_get_vara_double(ncid, varid, start, count, gps->pt[ids]));
694
695 NC(nc_inq_varid(ncid, "th", &varid));
696 NC(nc_get_vara_double(ncid, varid, start, count, &gps->th[ids]));
697 }
698
699 /* Close file... */
700 NC(nc_close(ncid));
701}

◆ read_met()

void read_met ( char *  filename,
met_t *  met 
)

Read meteorological data file.

Definition at line 705 of file libgps.c.

707 {
708
709 char tstr[10];
710
711 int dimid, ncid, varid, year, mon, day, hour;
712
713 size_t np, nx, ny;
714
715 /* Write info... */
716 printf("Read meteorological data: %s\n", filename);
717
718 /* Get time from filename... */
719 sprintf(tstr, "%.4s", &filename[strlen(filename) - 16]);
720 year = atoi(tstr);
721 sprintf(tstr, "%.2s", &filename[strlen(filename) - 11]);
722 mon = atoi(tstr);
723 sprintf(tstr, "%.2s", &filename[strlen(filename) - 8]);
724 day = atoi(tstr);
725 sprintf(tstr, "%.2s", &filename[strlen(filename) - 5]);
726 hour = atoi(tstr);
727 time2jsec(year, mon, day, hour, 0, 0, 0, &met->time);
728
729 /* Open netCDF file... */
730 NC(nc_open(filename, NC_NOWRITE, &ncid));
731
732 /* Get dimensions... */
733 NC(nc_inq_dimid(ncid, "lon", &dimid));
734 NC(nc_inq_dimlen(ncid, dimid, &nx));
735 if (nx > EX)
736 ERRMSG("Too many longitudes!");
737
738 NC(nc_inq_dimid(ncid, "lat", &dimid));
739 NC(nc_inq_dimlen(ncid, dimid, &ny));
740 if (ny > EY)
741 ERRMSG("Too many latitudes!");
742
743 NC(nc_inq_dimid(ncid, "lev", &dimid));
744 NC(nc_inq_dimlen(ncid, dimid, &np));
745 if (np > EP)
746 ERRMSG("Too many levels!");
747
748 /* Store dimensions... */
749 met->np = (int) np;
750 met->nx = (int) nx;
751 met->ny = (int) ny;
752
753 /* Get horizontal grid... */
754 NC(nc_inq_varid(ncid, "lon", &varid));
755 NC(nc_get_var_double(ncid, varid, met->lon));
756 NC(nc_inq_varid(ncid, "lat", &varid));
757 NC(nc_get_var_double(ncid, varid, met->lat));
758
759 /* Read meteorological data... */
760 read_met_help(ncid, "t", "T", met, met->t, 1.0);
761
762 /* Read pressure levels from file... */
763 NC(nc_inq_varid(ncid, "lev", &varid));
764 NC(nc_get_var_double(ncid, varid, met->p));
765 for (int ip = 0; ip < met->np; ip++)
766 met->p[ip] /= 100.;
767
768 /* Extrapolate data for lower boundary... */
769 read_met_extrapolate(met);
770
771 /* Check ordering of pressure levels... */
772 for (int ip = 1; ip < met->np; ip++)
773 if (met->p[ip - 1] < met->p[ip])
774 ERRMSG("Pressure levels must be descending!");
775
776 /* Create periodic boundary conditions... */
777 read_met_periodic(met);
778
779 /* Close file... */
780 NC(nc_close(ncid));
781}
read_met_extrapolate
void read_met_extrapolate(met_t *met)
Extrapolate meteorological data at lower boundary.
Definition: libgps.c:785
read_met_periodic
void read_met_periodic(met_t *met)
Create meteorological data with periodic boundary conditions.
Definition: libgps.c:838
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:806
EY
#define EY
Maximum number of latitudes for meteorological data.
Definition: libgps.h:97
EX
#define EX
Maximum number of longitudes for meteorological data.
Definition: libgps.h:94
EP
#define EP
Maximum number of pressure levels for meteorological data.
Definition: libgps.h:91
Here is the call graph for this function:

◆ read_met_extrapolate()

void read_met_extrapolate ( met_t *  met)

Extrapolate meteorological data at lower boundary.

Definition at line 785 of file libgps.c.

786 {
787
788 /* Loop over columns... */
789 for (int ix = 0; ix < met->nx; ix++)
790 for (int iy = 0; iy < met->ny; iy++) {
791
792 /* Find lowest valid data point... */
793 int ip0;
794 for (ip0 = met->np - 1; ip0 >= 0; ip0--)
795 if (!gsl_finite(met->t[ix][iy][ip0]))
796 break;
797
798 /* Extrapolate... */
799 for (int ip = ip0; ip >= 0; ip--)
800 met->t[ix][iy][ip] = met->t[ix][iy][ip + 1];
801 }
802}

◆ 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 at line 806 of file libgps.c.

812 {
813
814 static float help[EX * EY * EP];
815
816 int n = 0, varid;
817
818 /* Check if variable exists... */
819 if (nc_inq_varid(ncid, varname, &varid) != NC_NOERR)
820 if (nc_inq_varid(ncid, varname2, &varid) != NC_NOERR)
821 return;
822
823 /* Read data... */
824 NC(nc_get_var_float(ncid, varid, help));
825
826 /* Copy and check data... */
827 for (int ip = 0; ip < met->np; ip++)
828 for (int iy = 0; iy < met->ny; iy++)
829 for (int ix = 0; ix < met->nx; ix++) {
830 dest[ix][iy][ip] = scl * help[n++];
831 if (fabs(dest[ix][iy][ip] / scl) > 1e14)
832 dest[ix][iy][ip] = GSL_NAN;
833 }
834}

◆ read_met_periodic()

void read_met_periodic ( met_t *  met)

Create meteorological data with periodic boundary conditions.

Definition at line 838 of file libgps.c.

839 {
840
841 /* Check longitudes... */
842 if (!(fabs(met->lon[met->nx - 1] - met->lon[0]
843 + met->lon[1] - met->lon[0] - 360) < 0.01))
844 return;
845
846 /* Increase longitude counter... */
847 if ((++met->nx) > EX)
848 ERRMSG("Cannot create periodic boundary conditions!");
849
850 /* Set longitude... */
851 met->lon[met->nx - 1] = met->lon[met->nx - 2] + met->lon[1] - met->lon[0];
852
853 /* Loop over latitudes and pressure levels... */
854 for (int iy = 0; iy < met->ny; iy++)
855 for (int ip = 0; ip < met->np; ip++)
856 met->t[met->nx - 1][iy][ip] = met->t[0][iy][ip];
857}

◆ 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 at line 861 of file libgps.c.

868 {
869
870 /* Cubic spline interpolation... */
871 if (method == 1) {
872
873 /* Allocate... */
874 gsl_interp_accel *acc = gsl_interp_accel_alloc();
875 gsl_spline *s = gsl_spline_alloc(gsl_interp_cspline, (size_t) n);
876
877 /* Interpolate profile... */
878 gsl_spline_init(s, x, y, (size_t) n);
879 for (int i = 0; i < n2; i++)
880 if (x2[i] <= x[0])
881 y2[i] = y[0];
882 else if (x2[i] >= x[n - 1])
883 y2[i] = y[n - 1];
884 else
885 y2[i] = gsl_spline_eval(s, x2[i], acc);
886
887 /* Free... */
888 gsl_spline_free(s);
889 gsl_interp_accel_free(acc);
890 }
891
892 /* Linear interpolation... */
893 else {
894 for (int i = 0; i < n2; i++)
895 if (x2[i] <= x[0])
896 y2[i] = y[0];
897 else if (x2[i] >= x[n - 1])
898 y2[i] = y[n - 1];
899 else {
900 int idx = locate_irr(x, n, x2[i]);
901 y2[i] = LIN(x[idx], y[idx], x[idx + 1], y[idx + 1], x2[i]);
902 }
903 }
904}

◆ tropopause()

void tropopause ( gps_t *  gps)

Find tropopause height.

Definition at line 908 of file libgps.c.

909 {
910
911 /* Loop over profiles... */
912 for (int ids = 0; ids < gps->nds; ids++) {
913
914 /* Set default value... */
915 gps->th[ids] = GSL_NAN;
916
917 /* Set minimum altitude... */
918 const double zmin =
919 8 - 4 * fabs(cos((90 - gps->lat[ids][gps->nz[ids] / 2]) * M_PI / 180));
920
921 /* Search tropopause (WMO definition)... */
922 for (int iz = 0; iz < gps->nz[ids]; iz++)
923 if (gps->z[ids][iz] >= zmin && gps->z[ids][iz] <= 20.0) {
924 int okay = 1;
925 for (int iz2 = iz + 1; iz2 < gps->nz[ids]; iz2++)
926 if (gps->z[ids][iz2] - gps->z[ids][iz] <= 2.0)
927 if (!gsl_finite(gps->t[ids][iz]) ||
928 !gsl_finite(gps->t[ids][iz2]) ||
929 (gps->t[ids][iz2] - gps->t[ids][iz])
930 / (gps->z[ids][iz2] - gps->z[ids][iz]) < -2.0)
931 okay = 0;
932 if (okay) {
933 gps->th[ids] = gps->z[ids][iz];
934 break;
935 }
936 }
937 }
938}

◆ tropopause_spline()

void tropopause_spline ( gps_t *  gps,
int  met_tropo 
)

Find tropopause height using cubic spline interpolation.

Definition at line 942 of file libgps.c.

944 {
945
946 /* Loop over data sets... */
947#pragma omp parallel for default(shared)
948 for (int ids = 0; ids < gps->nds; ids++) {
949
950 /* Init... */
951 gps->tp[ids] = NAN;
952 gps->th[ids] = NAN;
953 gps->tt[ids] = NAN;
954 gps->tq[ids] = NAN;
955 gps->tlon[ids] = NAN;
956 gps->tlat[ids] = NAN;
957
958 /* Get vertical profiles... */
959 int nz = 0;
960 double h[NZ], t[NZ], z[NZ], q[NZ], lon[NZ], lat[NZ];
961 for (int iz = 0; iz < gps->nz[ids]; iz++)
962 if (gsl_finite(gps->p[ids][iz]) && gsl_finite(gps->t[ids][iz])
963 && gsl_finite(gps->z[ids][iz])
964 && gps->z[ids][iz] >= 4.0 && gps->z[ids][iz] <= 24.0) {
965 h[nz] = gps->z[ids][iz];
966 t[nz] = gps->t[ids][iz];
967 z[nz] = Z(gps->p[ids][iz]);
968 q[nz] = gps->wv[ids][iz];
969 lon[nz] = gps->lon[ids][iz];
970 lat[nz] = gps->lat[ids][iz];
971 if (nz > 0 && z[nz] <= z[nz - 1])
972 ERRMSG("Profiles must be ascending!");
973 if ((++nz) >= NZ)
974 ERRMSG("Too many height levels!");
975 }
976 if (z[0] > 4.5 || z[nz - 1] < 23.5)
977 WARN("Vertical profile is incomplete!");
978
979 /* Set grid for spline interpolation... */
980 double h2[200], p2[200], t2[200], z2[200], q2[200];
981 for (int iz = 0; iz <= 190; iz++) {
982 z2[iz] = 4.5 + 0.1 * iz;
983 p2[iz] = P(z2[iz]);
984 }
985
986 /* Interpolate temperature and geopotential height profiles... */
987 spline(z, t, nz, z2, t2, 191, 1);
988 spline(z, h, nz, z2, h2, 191, 1);
989 spline(z, q, nz, z2, q2, 191, 1);
990
991 /* Use cold point... */
992 if (met_tropo == 2) {
993
994 /* Find minimum... */
995 int iz = (int) gsl_stats_min_index(t2, 1, 171);
996 if (iz > 0 && iz < 170) {
997 gps->tp[ids] = p2[iz];
998 gps->th[ids] = h2[iz];
999 gps->tt[ids] = t2[iz];
1000 gps->tq[ids] = q2[iz];
1001 }
1002 }
1003
1004 /* Use WMO definition... */
1005 else if (met_tropo == 3 || met_tropo == 4) {
1006
1007 /* Find 1st tropopause... */
1008 int iz;
1009 for (iz = 0; iz <= 170; iz++) {
1010 int found = 1;
1011 for (int iz2 = iz + 1; iz2 <= iz + 20; iz2++)
1012 if (LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]) > 2.0) {
1013 found = 0;
1014 break;
1015 }
1016 if (found) {
1017 if (iz > 0 && iz < 170) {
1018 gps->tp[ids] = p2[iz];
1019 gps->th[ids] = h2[iz];
1020 gps->tt[ids] = t2[iz];
1021 gps->tq[ids] = q2[iz];
1022 }
1023 break;
1024 }
1025 }
1026
1027 /* Find 2nd tropopause... */
1028 if (met_tropo == 4) {
1029
1030 /* Init... */
1031 gps->tp[ids] = NAN;
1032 gps->th[ids] = NAN;
1033 gps->tt[ids] = NAN;
1034 gps->tq[ids] = NAN;
1035
1036 /* Check layers... */
1037 for (; iz <= 170; iz++) {
1038 int found = 1;
1039 for (int iz2 = iz + 1; iz2 <= iz + 10; iz2++)
1040 if (LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]) < 3.0) {
1041 found = 0;
1042 break;
1043 }
1044 if (found)
1045 break;
1046 }
1047 for (; iz <= 170; iz++) {
1048 int found = 1;
1049 for (int iz2 = iz + 1; iz2 <= iz + 20; iz2++)
1050 if (LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]) > 2.0) {
1051 found = 0;
1052 break;
1053 }
1054 if (found) {
1055 if (iz > 0 && iz < 170) {
1056 gps->tp[ids] = p2[iz];
1057 gps->th[ids] = h2[iz];
1058 gps->tt[ids] = t2[iz];
1059 gps->tq[ids] = q2[iz];
1060 }
1061 break;
1062 }
1063 }
1064 }
1065 }
1066
1067 /* Find tropopause longitude and latitude... */
1068 if (gsl_finite(gps->th[ids]))
1069 for (int iz = 0; iz < nz - 1; iz++)
1070 if (gps->th[ids] >= h[iz] && gps->th[ids] < h[iz + 1]) {
1071 gps->tlon[ids] = lon[iz];
1072 gps->tlat[ids] = lat[iz];
1073 break;
1074 }
1075 }
1076}
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:861
LAPSE
#define LAPSE(p1, t1, p2, t2)
Calculate lapse rate.
Definition: libgps.h:129
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
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::tp
double tp[NDS]
Tropopause pressure [hPa].
Definition: libgps.h:189
gps_t::tlat
double tlat[NDS]
Tropopause latitude [deg].
Definition: libgps.h:201
gps_t::tq
double tq[NDS]
Tropopause water vapor [ppmv].
Definition: libgps.h:195
Here is the call graph for this function:

◆ write_gps()

void write_gps ( char *  filename,
gps_t *  gps 
)

Write GPS-RO data file.

Definition at line 1080 of file libgps.c.

1082 {
1083
1084 static double help[NDS * NZ];
1085
1086 int ncid, dimid[2], time_id, z_id, lon_id, lat_id, p_id, t_id,
1087 pt_id, wv_id, th_id, nzmax = 0;
1088
1089 /* Create netCDF file... */
1090 printf("Write GPS-RO file: %s\n", filename);
1091 NC(nc_create(filename, NC_CLOBBER, &ncid));
1092
1093 /* Set dimensions... */
1094 NC(nc_def_dim(ncid, "NDS", (size_t) gps->nds, &dimid[0]));
1095 for (int ids = 0; ids < gps->nds; ids++)
1096 nzmax = GSL_MAX(nzmax, gps->nz[ids]);
1097 NC(nc_def_dim(ncid, "NZ", (size_t) nzmax, &dimid[1]));
1098
1099 /* Add variables... */
1100 add_var(ncid, "time", "s", "time (seconds since 2000-01-01T00:00Z)",
1101 NC_DOUBLE, dimid, &time_id, 1);
1102 add_var(ncid, "z", "km", "altitude", NC_FLOAT, dimid, &z_id, 2);
1103 add_var(ncid, "lon", "deg", "longitude", NC_FLOAT, dimid, &lon_id, 2);
1104 add_var(ncid, "lat", "deg", "latitude", NC_FLOAT, dimid, &lat_id, 2);
1105 add_var(ncid, "p", "hPa", "pressure", NC_FLOAT, dimid, &p_id, 2);
1106 add_var(ncid, "t", "K", "temperature", NC_FLOAT, dimid, &t_id, 2);
1107 add_var(ncid, "wv", "ppv", "water vapor volume mixing ratio",
1108 NC_FLOAT, dimid, &wv_id, 2);
1109 add_var(ncid, "pt", "K", "temperature perturbation",
1110 NC_FLOAT, dimid, &pt_id, 2);
1111 add_var(ncid, "th", "km", "tropopause height", NC_FLOAT, dimid, &th_id, 1);
1112
1113 /* Leave define mode... */
1114 NC(nc_enddef(ncid));
1115
1116 /* Write data... */
1117 NC(nc_put_var_double(ncid, time_id, gps->time));
1118 NC(nc_put_var_double(ncid, th_id, gps->th));
1119 for (int ids = 0; ids < gps->nds; ids++)
1120 for (int iz = 0; iz < gps->nz[ids]; iz++)
1121 help[ids * nzmax + iz] = gps->z[ids][iz];
1122 NC(nc_put_var_double(ncid, z_id, help));
1123 for (int ids = 0; ids < gps->nds; ids++)
1124 for (int iz = 0; iz < gps->nz[ids]; iz++)
1125 help[ids * nzmax + iz] = gps->lon[ids][iz];
1126 NC(nc_put_var_double(ncid, lon_id, help));
1127 for (int ids = 0; ids < gps->nds; ids++)
1128 for (int iz = 0; iz < gps->nz[ids]; iz++)
1129 help[ids * nzmax + iz] = gps->lat[ids][iz];
1130 NC(nc_put_var_double(ncid, lat_id, help));
1131 for (int ids = 0; ids < gps->nds; ids++)
1132 for (int iz = 0; iz < gps->nz[ids]; iz++)
1133 help[ids * nzmax + iz] = gps->p[ids][iz];
1134 NC(nc_put_var_double(ncid, p_id, help));
1135 for (int ids = 0; ids < gps->nds; ids++)
1136 for (int iz = 0; iz < gps->nz[ids]; iz++)
1137 help[ids * nzmax + iz] = gps->t[ids][iz];
1138 NC(nc_put_var_double(ncid, t_id, help));
1139 for (int ids = 0; ids < gps->nds; ids++)
1140 for (int iz = 0; iz < gps->nz[ids]; iz++)
1141 help[ids * nzmax + iz] = gps->wv[ids][iz];
1142 NC(nc_put_var_double(ncid, wv_id, help));
1143 for (int ids = 0; ids < gps->nds; ids++)
1144 for (int iz = 0; iz < gps->nz[ids]; iz++)
1145 help[ids * nzmax + iz] = gps->pt[ids][iz];
1146 NC(nc_put_var_double(ncid, pt_id, help));
1147
1148 /* Close file... */
1149 NC(nc_close(ncid));
1150}
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
Here is the call graph for this function:
Generated by doxygen 1.9.4