GPS Code Collection
Macros | Functions
variance.c File Reference

Estimate gravity wave variances. More...

#include "libgps.h"

Go to the source code of this file.

Macros

#define GX   360
 
#define GY   180
 
#define GZ   50
 

Functions

int main (int argc, char *argv[])
 

Detailed Description

Estimate gravity wave variances.

Definition in file variance.c.

Macro Definition Documentation

◆ GX

#define GX   360

Definition at line 33 of file variance.c.

◆ GY

#define GY   180

Definition at line 36 of file variance.c.

◆ GZ

#define GZ   50

Definition at line 39 of file variance.c.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 45 of file variance.c.

47 {
48
49 gps_t *gps;
50
51 FILE *out;
52
53 static double lon0, lon1, lat0, lat1, z0, z1, mean[GX][GY][GZ],
54 min[GX][GY][GZ], max[GX][GY][GZ], var[GX][GY][GZ],
55 mtime[GX][GY], glon[GX], glat[GY], gz[GZ], thmean[GX][GY],
56 tmean[GX][GY][GZ], twmean[GX][GY], se[NZ], sz[NZ], tw, w, wmax, wsum;
57
58 static int iarg, ids, idx, ix, iy, iz, iz2,
59 nx, ny, nz, n[GX][GY][GZ], np[GX][GY], sn;
60
61 /* Allocate... */
62 ALLOC(gps, gps_t, 1);
63
64 /* Check arguments... */
65 if (argc < 5)
66 ERRMSG("Give parameters: <ctl> <var.tab> <sens.tab> "
67 "<gps1.nc> [<gps2.nx> ...]");
68
69 /* Get control parameters... */
70 z0 = scan_ctl(argc, argv, "Z0", -1, "0", NULL);
71 z1 = scan_ctl(argc, argv, "Z1", -1, "60", NULL);
72 nz = (int) scan_ctl(argc, argv, "NZ", -1, "6", NULL);
73 lon0 = scan_ctl(argc, argv, "LON0", -1, "-180", NULL);
74 lon1 = scan_ctl(argc, argv, "LON1", -1, "180", NULL);
75 nx = (int) scan_ctl(argc, argv, "NX", -1, "72", NULL);
76 lat0 = scan_ctl(argc, argv, "LAT0", -1, "-90", NULL);
77 lat1 = scan_ctl(argc, argv, "LAT1", -1, "90", NULL);
78 ny = (int) scan_ctl(argc, argv, "NY", -1, "36", NULL);
79
80 /* Check grid dimensions... */
81 if (nx < 1 || nx > GX)
82 ERRMSG("Set 1 <= GX <= MAX!");
83 if (ny < 1 || ny > GY)
84 ERRMSG("Set 1 <= GY <= MAX!");
85 if (nz < 1 || nz > GZ)
86 ERRMSG("Set 1 <= GZ <= MAX!");
87
88 /* Read vertical sensitivity function... */
89 if (argv[3][0] != '-') {
90 read_shape(argv[3], sz, se, &sn);
91 if (sn > NZ)
92 ERRMSG("Too many data points!");
93 }
94
95 /* Loop over data files... */
96 for (iarg = 4; iarg < argc; iarg++) {
97
98 /* Read gps data... */
99 FILE *in;
100 if (!(in = fopen(argv[iarg], "r")))
101 continue;
102 else {
103 fclose(in);
104 read_gps(argv[iarg], gps);
105 }
106
107 /* Loop over profiles... */
108 for (ids = 0; ids < gps->nds; ids++) {
109
110 /* Check tropopause height... */
111 if (!gsl_finite(gps->th[ids]))
112 continue;
113
114 /* Multiply with vertical sensitivity function... */
115 if (argv[3][0] != '-') {
116 tw = wsum = wmax = 0;
117 for (iz2 = 0; iz2 < gps->nz[ids]; iz2++) {
118 if (gps->z[ids][iz2] < sz[0] || gps->z[ids][iz2] > sz[sn - 1])
119 w = 0;
120 else {
121 idx = locate_irr(sz, sn, gps->z[ids][iz2]);
122 w =
123 LIN(sz[idx], se[idx], sz[idx + 1], se[idx + 1],
124 gps->z[ids][iz2]);
125 }
126 if (gsl_finite(gps->t[ids][iz2]) && gps->pt[ids][iz2]) {
127 tw += w * gps->t[ids][iz2];
128 wsum += w;
129 gps->pt[ids][iz2] *= w;
130 wmax = GSL_MAX(w, wmax);
131 }
132 }
133 tw /= wsum;
134 for (iz2 = 0; iz2 < gps->nz[ids]; iz2++)
135 gps->pt[ids][iz2] /= wmax;
136 }
137
138 /* Get grid indices... */
139 ix = (int) ((gps->lon[ids][gps->nz[ids] / 2] - lon0)
140 / (lon1 - lon0) * (double) nx);
141 iy = (int) ((gps->lat[ids][gps->nz[ids] / 2] - lat0)
142 / (lat1 - lat0) * (double) ny);
143 if (ix < 0 || ix >= nx || iy < 0 || iy >= ny)
144 continue;
145
146 /* Get mean time and tropopause height... */
147 mtime[ix][iy] += gps->time[ids];
148 thmean[ix][iy] += gps->th[ids];
149 twmean[ix][iy] += tw;
150 np[ix][iy]++;
151
152 /* Loop over altitudes... */
153 for (iz2 = 0; iz2 < gps->nz[ids]; iz2++) {
154
155 /* Get grid indices... */
156 iz = (int) ((gps->z[ids][iz2] - z0)
157 / (z1 - z0) * (double) nz);
158 if (iz < 0 || iz >= nz)
159 continue;
160
161 /* Check data... */
162 if (!gsl_finite(gps->t[ids][iz2]) || !gsl_finite(gps->pt[ids][iz2]))
163 continue;
164
165 /* Get statistics of perturbations... */
166 tmean[ix][iy][iz] += gps->t[ids][iz2];
167 mean[ix][iy][iz] += gps->pt[ids][iz2];
168 var[ix][iy][iz] += gsl_pow_2(gps->pt[ids][iz2]);
169 max[ix][iy][iz] = GSL_MAX(max[ix][iy][iz], gps->pt[ids][iz2]);
170 min[ix][iy][iz] = GSL_MIN(min[ix][iy][iz], gps->pt[ids][iz2]);
171 n[ix][iy][iz]++;
172 }
173 }
174 }
175
176 /* Analyze results... */
177 for (ix = 0; ix < nx; ix++)
178 for (iy = 0; iy < ny; iy++) {
179
180 /* Get mean time and tropopause height... */
181 if (np[ix][iy] > 0) {
182 mtime[ix][iy] /= (double) np[ix][iy];
183 thmean[ix][iy] /= (double) np[ix][iy];
184 twmean[ix][iy] /= (double) np[ix][iy];
185 } else {
186 mtime[ix][iy] = GSL_NAN;
187 thmean[ix][iy] = GSL_NAN;
188 twmean[ix][iy] = GSL_NAN;
189 }
190
191 /* Loop over altitudes... */
192 for (iz = 0; iz < nz; iz++) {
193
194 /* Get geolocation... */
195 gz[iz] = z0 + (iz + 0.5) / (double) nz *(
196 z1 - z0);
197 glon[ix] = lon0 + (ix + 0.5) / (double) nx *(
198 lon1 - lon0);
199 glat[iy] = lat0 + (iy + 0.5) / (double) ny *(
200 lat1 - lat0);
201
202 /* Get mean perturbation and variance... */
203 if (n[ix][iy][iz] > 0) {
204 tmean[ix][iy][iz]
205 /= (double) n[ix][iy][iz];
206 mean[ix][iy][iz]
207 /= (double) n[ix][iy][iz];
208 var[ix][iy][iz]
209 = var[ix][iy][iz] / (double) n[ix][iy][iz]
210 - gsl_pow_2(mean[ix][iy][iz]);
211 } else {
212 tmean[ix][iy][iz] = GSL_NAN;
213 mean[ix][iy][iz] = GSL_NAN;
214 var[ix][iy][iz] = GSL_NAN;
215 min[ix][iy][iz] = GSL_NAN;
216 max[ix][iy][iz] = GSL_NAN;
217 }
218 }
219 }
220
221 /* Create file... */
222 printf("Write variance statistics: %s\n", argv[2]);
223 if (!(out = fopen(argv[2], "w")))
224 ERRMSG("Cannot create file!");
225
226 /* Write header... */
227 fprintf(out,
228 "# $1 = time [s]\n"
229 "# $2 = altitude [km]\n"
230 "# $3 = longitude [deg]\n"
231 "# $4 = latitude [deg]\n"
232 "# $5 = number of profiles\n"
233 "# $6 = number of data points\n"
234 "# $7 = mean perturbation [K]\n"
235 "# $8 = minimum perturbation [K]\n"
236 "# $9 = maximum perturbation [K]\n"
237 "# $10 = variance [K^2]\n"
238 "# $11 = mean temperature [K]\n"
239 "# $12 = mean weighted temperature [K]\n"
240 "# $13 = mean tropopause height [km]\n");
241
242 /* Write results... */
243 for (iz = 0; iz < nz; iz++) {
244 for (iy = 0; iy < ny; iy++) {
245 if (iy == 0 || nx > 1)
246 fprintf(out, "\n");
247 for (ix = 0; ix < nx; ix++)
248 fprintf(out, "%.2f %g %g %g %d %d %g %g %g %g %g %g %g\n",
249 mtime[ix][iy], gz[iz], glon[ix], glat[iy],
250 np[ix][iy], n[ix][iy][iz], mean[ix][iy][iz],
251 min[ix][iy][iz], max[ix][iy][iz], var[ix][iy][iz],
252 tmean[ix][iy][iz], twmean[ix][iy], thmean[ix][iy]);
253 }
254 }
255
256 /* Close file... */
257 fclose(out);
258
259 /* Free... */
260 free(gps);
261
262 return EXIT_SUCCESS;
263}
read_gps
void read_gps(char *filename, gps_t *gps)
Read GPS-RO data file.
Definition: libgps.c:632
NZ
#define NZ
Maximum number of altitudes per GPS-RO profile.
Definition: libgps.h:103
gps_t
GPS-RO profile data.
Definition: libgps.h:153
gps_t::time
double time[NDS]
Time (seconds since 2000-01-01T00:00Z).
Definition: libgps.h:162
gps_t::pt
double pt[NDS][NZ]
Temperature perturbation [K].
Definition: libgps.h:183
gps_t::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::th
double th[NDS]
Tropopause height [km].
Definition: libgps.h:186
gps_t::z
double z[NDS][NZ]
Altitude [km].
Definition: libgps.h:165
gps_t::lat
double lat[NDS][NZ]
Latitude [deg].
Definition: libgps.h:171
GY
#define GY
Definition: variance.c:36
GX
#define GX
Definition: variance.c:33
GZ
#define GZ
Definition: variance.c:39
Here is the call graph for this function:
Generated on Tue Aug 19 2025 18:04:34 for GPS Code Collection by doxygen 1.9.4