GPS Code Collection
Data Structures | Macros | Functions
libgps.h File Reference

GPS Code Collection library declarations. More...

#include <netcdf.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_spline.h>
#include "jurassic.h"

Go to the source code of this file.

Data Structures

struct  gps_t
 GPS-RO profile data. More...
 
struct  met_t
 Meteorological data. More...
 

Macros

#define EP   73
 Maximum number of pressure levels for meteorological data. More...
 
#define EX   721
 Maximum number of longitudes for meteorological data. More...
 
#define EY   361
 Maximum number of latitudes for meteorological data. More...
 
#define NDS   10000
 Maximum number of GPS-RO profiles. More...
 
#define NZ   5000
 Maximum number of altitudes per GPS-RO profile. More...
 
#define MA   28.9644
 Molar mass of dry air [g/mol]. More...
 
#define RA   (1e3 * RI / MA)
 Specific gas constant of dry air [J/(kg K)]. More...
 
#define RI   8.3144598
 Ideal gas constant [J/(mol K)]. More...
 
#define LAPSE(p1, t1, p2, t2)
 Calculate lapse rate. More...
 
#define NC(cmd)
 Execute netCDF library command and check result. More...
 
#define P(z)    (P0 * exp(-(z) / H0))
 Compute pressure at given altitude. More...
 
#define Z(p)    (H0 * log(P0 / (p)))
 Convert pressure to altitude. More...
 

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 declarations.

Definition in file libgps.h.

Macro Definition Documentation

◆ EP

#define EP   73

Maximum number of pressure levels for meteorological data.

Definition at line 91 of file libgps.h.

◆ EX

#define EX   721

Maximum number of longitudes for meteorological data.

Definition at line 94 of file libgps.h.

◆ EY

#define EY   361

Maximum number of latitudes for meteorological data.

Definition at line 97 of file libgps.h.

◆ NDS

#define NDS   10000

Maximum number of GPS-RO profiles.

Definition at line 100 of file libgps.h.

◆ NZ

#define NZ   5000

Maximum number of altitudes per GPS-RO profile.

Definition at line 103 of file libgps.h.

◆ MA

#define MA   28.9644

Molar mass of dry air [g/mol].

Definition at line 111 of file libgps.h.

◆ RA

#define RA   (1e3 * RI / MA)

Specific gas constant of dry air [J/(kg K)].

Definition at line 116 of file libgps.h.

◆ RI

#define RI   8.3144598

Ideal gas constant [J/(mol K)].

Definition at line 121 of file libgps.h.

◆ LAPSE

#define LAPSE (   p1,
  t1,
  p2,
  t2 
)
Value:
(1e3 * G0 / RA * ((t2) - (t1)) / ((t2) + (t1)) \
* ((p2) + (p1)) / ((p2) - (p1)))
RA
#define RA
Specific gas constant of dry air [J/(kg K)].
Definition: libgps.h:116

Calculate lapse rate.

Definition at line 129 of file libgps.h.

◆ NC

#define NC (   cmd)
Value:
{ \
int nc_result=(cmd); \
if(nc_result!=NC_NOERR) \
ERRMSG("%s", nc_strerror(nc_result)); \
}

Execute netCDF library command and check result.

Definition at line 134 of file libgps.h.

◆ P

#define P (   z)     (P0 * exp(-(z) / H0))

Compute pressure at given altitude.

Definition at line 141 of file libgps.h.

◆ Z

#define Z (   p)     (H0 * log(P0 / (p)))

Convert pressure to altitude.

Definition at line 145 of file libgps.h.

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 double t;
65
66 int ids, iz;
67
68 /* Allocate... */
69 ALLOC(met0, met_t, 1);
70 ALLOC(met1, met_t, 1);
71
72 /* Loop over profiles... */
73 for (ids = 0; ids < gps->nds; ids++) {
74
75 /* Loop over altitudes... */
76 for (iz = 0; iz < gps->nz[ids]; iz++) {
77
78 /* Get meteorological data... */
79 get_met(metbase, dt_met, gps->time[ids], met0, met1);
80
81 /* Interpolate meteorological data... */
82 intpol_met_time(met0, met1, gps->time[ids], gps->p[ids][iz],
83 gps->lon[ids][iz], gps->lat[ids][iz], &t);
84
85 /* Set perturbation... */
86 gps->pt[ids][iz] = gps->t[ids][iz] - t;
87 }
88 }
89
90 /* Free... */
91 free(met0);
92 free(met1);
93}
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:316
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:197
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 97 of file libgps.c.

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

◆ grid_gps()

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

Interpolate GPS data to regular altitude grid.

Definition at line 140 of file libgps.c.

144 {
145
146 double lat[NZ], lon[NZ], p[NZ], pt[NZ], t[NZ], wv[NZ], z[NZ];
147
148 int ids, iz, iz2;
149
150 /* Check number of altitudes... */
151 if (nz > NZ)
152 ERRMSG("Too many altitudes!");
153
154 /* Loop over profiles... */
155 for (ids = 0; ids < gps->nds; ids++) {
156
157 /* Loop over altitudes... */
158 for (iz = 0; iz < nz; iz++) {
159
160 /* Set altitude... */
161 z[iz] = LIN(0.0, zmin, nz - 1.0, zmax, (double) iz);
162
163 /* Get index... */
164 iz2 = locate_irr(gps->z[ids], gps->nz[ids], z[iz]);
165
166 /* Interpolate... */
167 lon[iz] = LIN(gps->z[ids][iz2], gps->lon[ids][iz2],
168 gps->z[ids][iz2 + 1], gps->lon[ids][iz2 + 1], z[iz]);
169 lat[iz] = LIN(gps->z[ids][iz2], gps->lat[ids][iz2],
170 gps->z[ids][iz2 + 1], gps->lat[ids][iz2 + 1], z[iz]);
171 p[iz] = LIN(gps->z[ids][iz2], gps->p[ids][iz2],
172 gps->z[ids][iz2 + 1], gps->p[ids][iz2 + 1], z[iz]);
173 t[iz] = LIN(gps->z[ids][iz2], gps->t[ids][iz2],
174 gps->z[ids][iz2 + 1], gps->t[ids][iz2 + 1], z[iz]);
175 wv[iz] = LIN(gps->z[ids][iz2], gps->wv[ids][iz2],
176 gps->z[ids][iz2 + 1], gps->wv[ids][iz2 + 1], z[iz]);
177 pt[iz] = LIN(gps->z[ids][iz2], gps->pt[ids][iz2],
178 gps->z[ids][iz2 + 1], gps->pt[ids][iz2 + 1], z[iz]);
179 }
180
181 /* Copy data... */
182 gps->nz[ids] = nz;
183 for (iz = 0; iz < nz; iz++) {
184 gps->z[ids][iz] = z[iz];
185 gps->lon[ids][iz] = lon[iz];
186 gps->lat[ids][iz] = lat[iz];
187 gps->p[ids][iz] = p[iz];
188 gps->t[ids][iz] = t[iz];
189 gps->wv[ids][iz] = wv[iz];
190 gps->pt[ids][iz] = pt[iz];
191 }
192 }
193}
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 197 of file libgps.c.

202 {
203
204 char filename[LEN];
205
206 static int init;
207
208 /* Init... */
209 if (!init) {
210 init = 1;
211
212 get_met_help(t, -1, metbase, dt_met, filename);
213 read_met(filename, met0);
214
215 get_met_help(t + 1.0, 1, metbase, dt_met, filename);
216 read_met(filename, met1);
217 }
218
219 /* Read new data... */
220 if (t > met1->time) {
221 memcpy(met0, met1, sizeof(met_t));
222 get_met_help(t, 1, metbase, dt_met, filename);
223 read_met(filename, met1);
224 }
225}
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:229
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 229 of file libgps.c.

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

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

263 {
264
265 double aux00, aux01, aux10, aux11;
266
267 /* Interpolate vertically... */
268 aux00 = wp * (array[ix][iy][ip] - array[ix][iy][ip + 1])
269 + array[ix][iy][ip + 1];
270 aux01 = wp * (array[ix][iy + 1][ip] - array[ix][iy + 1][ip + 1])
271 + array[ix][iy + 1][ip + 1];
272 aux10 = wp * (array[ix + 1][iy][ip] - array[ix + 1][iy][ip + 1])
273 + array[ix + 1][iy][ip + 1];
274 aux11 = wp * (array[ix + 1][iy + 1][ip] - array[ix + 1][iy + 1][ip + 1])
275 + array[ix + 1][iy + 1][ip + 1];
276
277 /* Interpolate horizontally... */
278 aux00 = wy * (aux00 - aux01) + aux01;
279 aux11 = wy * (aux10 - aux11) + aux11;
280 *var = wx * (aux00 - aux11) + aux11;
281}

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

290 {
291
292 double wp, wx, wy;
293
294 int ip, ix, iy;
295
296 /* Check longitude... */
297 if (met->lon[met->nx - 1] > 180 && lon < 0)
298 lon += 360;
299
300 /* Get indices... */
301 ip = locate_irr(met->p, met->np, p);
302 ix = locate_reg(met->lon, met->nx, lon);
303 iy = locate_reg(met->lat, met->ny, lat);
304
305 /* Get weights... */
306 wp = (met->p[ip + 1] - p) / (met->p[ip + 1] - met->p[ip]);
307 wx = (met->lon[ix + 1] - lon) / (met->lon[ix + 1] - met->lon[ix]);
308 wy = (met->lat[iy + 1] - lat) / (met->lat[iy + 1] - met->lat[iy]);
309
310 /* Interpolate... */
311 intpol_met_3d(met->t, ip, ix, iy, wp, wx, wy, t);
312}
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:255
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 316 of file libgps.c.

323 {
324
325 double t0, t1, wt;
326
327 /* Spatial interpolation... */
328 intpol_met_space(met0, p, lon, lat, &t0);
329 intpol_met_space(met1, p, lon, lat, &t1);
330
331 /* Get weighting factor... */
332 wt = (met1->time - ts) / (met1->time - met0->time);
333
334 /* Interpolate... */
335 *t = wt * (t0 - t1) + t1;
336}
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:285
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 340 of file libgps.c.

342 {
343
344 double ham[NZ], wsum;
345
346 int ids, iham, iz, nham;
347
348 /* Loop over profiles... */
349 for (ids = 0; ids < gps->nds; ids++) {
350
351 /* Calculate Hamming window coefficients... */
352 nham = (int) (dz / fabs((gps->z[ids][0] - gps->z[ids][gps->nz[ids] - 1])
353 / (gps->nz[ids] - 1.0)) + 0.5);
354 nham = GSL_MAX(GSL_MIN(nham, NZ), 2);
355 for (iham = 0; iham < nham; iham++)
356 ham[iham] = 0.54 + 0.46 * cos(M_PI * iham / (nham - 1.0));
357
358 /* Loop over altitudes... */
359 for (iz = 0; iz < gps->nz[ids]; iz++) {
360
361 /* Initialize... */
362 gps->pt[ids][iz] = ham[0] * gps->t[ids][iz];
363 wsum = ham[0];
364
365 /* Loop over filter window... */
366 for (iham = 1; iham < nham; iham++) {
367
368 /* Check array range... */
369 if (iz - iham < 0 || iz + iham >= gps->nz[ids])
370 continue;
371
372 /* Check temperature value... */
373 if (!gsl_finite(gps->t[ids][iz - iham]) ||
374 !gsl_finite(gps->t[ids][iz + iham]))
375 continue;
376
377 /* Check for tropopause... */
378 if (gsl_finite(gps->th[ids]) && gps->th[ids] > 0)
379 if ((gps->z[ids][iz] >= gps->th[ids]
380 && gps->z[ids][iz - iham] < gps->th[ids])
381 || (gps->z[ids][iz] <= gps->th[ids]
382 && gps->z[ids][iz + iham] > gps->th[ids]))
383 continue;
384
385 /* Apply Hamming filter... */
386 gps->pt[ids][iz]
387 += ham[iham] * (gps->t[ids][iz - iham] + gps->t[ids][iz + iham]);
388 wsum += 2 * ham[iham];
389 }
390
391 /* Calculate perturbation... */
392 gps->pt[ids][iz] = gps->t[ids][iz] - gps->pt[ids][iz] / wsum;
393 }
394 }
395}
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 399 of file libgps.c.

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

◆ poly()

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

Remove polynomial fit from perturbation profile.

Definition at line 454 of file libgps.c.

458 {
459
460 double bg[NZ];
461
462 int ids, iz;
463
464 /* Loop over profiles... */
465 for (ids = 0; ids < gps->nds; ids++) {
466
467 /* Set profile... */
468 for (iz = 0; iz < gps->nz[ids]; iz++)
469 bg[iz] = gps->pt[ids][iz];
470
471 /* Polynomial interpolation... */
472 poly_help(gps->z[ids], bg, gps->nz[ids], dim, zmin, zmax);
473
474 /* Remove background... */
475 for (iz = 0; iz < gps->nz[ids]; iz++)
476 gps->pt[ids][iz] -= bg[iz];
477 }
478}
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:482
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 482 of file libgps.c.

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

◆ read_gps_prof()

void read_gps_prof ( char *  filename,
gps_t *  gps 
)

Read GPS-RO profile.

Definition at line 541 of file libgps.c.

543 {
544
545 char bad[10];
546
547 double t0, t1, zmin = 1e100, zmax = -1e100;
548
549 int ncid, dimid, varid;
550
551 size_t iz, nz;
552
553 /* Open netCDF file... */
554 printf("Read GPS-RO profile: %s\n", filename);
555 NC(nc_open(filename, NC_NOWRITE, &ncid));
556
557 /* Get dimensions... */
558 NC(nc_inq_dimid(ncid, "MSL_alt", &dimid));
559 NC(nc_inq_dimlen(ncid, dimid, &nz));
560 gps->nz[gps->nds] = (int) nz;
561 if (nz > NZ)
562 ERRMSG("Too many altitudes!");
563
564 /* Check data quality flag... */
565 NC(nc_get_att_text(ncid, NC_GLOBAL, "bad", bad));
566 if (bad[0] != '0') {
567 NC(nc_close(ncid));
568 return;
569 }
570
571 /* Get time... */
572 NC(nc_get_att_double(ncid, NC_GLOBAL, "start_time", &t0));
573 NC(nc_get_att_double(ncid, NC_GLOBAL, "stop_time", &t1));
574 gps->time[gps->nds] = 0.5 * (t0 + t1) - 630720000.0;
575
576 /* Get data... */
577 NC(nc_inq_varid(ncid, "MSL_alt", &varid));
578 NC(nc_get_var_double(ncid, varid, gps->z[gps->nds]));
579 NC(nc_inq_varid(ncid, "Lon", &varid));
580 NC(nc_get_var_double(ncid, varid, gps->lon[gps->nds]));
581 NC(nc_inq_varid(ncid, "Lat", &varid));
582 NC(nc_get_var_double(ncid, varid, gps->lat[gps->nds]));
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 (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 (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 (iz = 0; iz < nz; iz++)
616 gps->t[gps->nds][iz] += 273.15;
617
618 /* Convert water vapor... */
619 for (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 ip, 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 (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 (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:837
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:805
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 int ip, ip0, ix, iy;
787
788 /* Loop over columns... */
789 for (ix = 0; ix < met->nx; ix++)
790 for (iy = 0; iy < met->ny; iy++) {
791
792 /* Find lowest valid data point... */
793 for (ip0 = met->np - 1; ip0 >= 0; ip0--)
794 if (!gsl_finite(met->t[ix][iy][ip0]))
795 break;
796
797 /* Extrapolate... */
798 for (ip = ip0; ip >= 0; ip--)
799 met->t[ix][iy][ip] = met->t[ix][iy][ip + 1];
800 }
801}

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

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

◆ read_met_periodic()

void read_met_periodic ( met_t *  met)

Create meteorological data with periodic boundary conditions.

Definition at line 837 of file libgps.c.

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

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

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

◆ tropopause()

void tropopause ( gps_t *  gps)

Find tropopause height.

Definition at line 909 of file libgps.c.

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

◆ tropopause_spline()

void tropopause_spline ( gps_t *  gps,
int  met_tropo 
)

Find tropopause height using cubic spline interpolation.

Definition at line 947 of file libgps.c.

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

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