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)
 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:703
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 
)

Read GPS-RO profile.

Definition at line 528 of file libgps.c.

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

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

◆ read_met()

void read_met ( char *  filename,
met_t *  met 
)

Read meteorological data file.

Definition at line 703 of file libgps.c.

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

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

◆ 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 804 of file libgps.c.

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

◆ read_met_periodic()

void read_met_periodic ( met_t *  met)

Create meteorological data with periodic boundary conditions.

Definition at line 836 of file libgps.c.

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

◆ 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 859 of file libgps.c.

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

◆ tropopause()

void tropopause ( gps_t *  gps)

Find tropopause height.

Definition at line 906 of file libgps.c.

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

◆ tropopause_spline()

void tropopause_spline ( gps_t *  gps,
int  met_tropo 
)

Find tropopause height using cubic spline interpolation.

Definition at line 940 of file libgps.c.

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

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