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 *in, *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 if (!(in = fopen(argv[iarg], "r")))
100 continue;
101 else {
102 fclose(in);
103 read_gps(argv[iarg], gps);
104 }
105
106 /* Loop over profiles... */
107 for (ids = 0; ids < gps->nds; ids++) {
108
109 /* Check tropopause height... */
110 if (!gsl_finite(gps->th[ids]))
111 continue;
112
113 /* Multiply with vertical sensitivity function... */
114 if (argv[3][0] != '-') {
115 tw = wsum = wmax = 0;
116 for (iz2 = 0; iz2 < gps->nz[ids]; iz2++) {
117 if (gps->z[ids][iz2] < sz[0] || gps->z[ids][iz2] > sz[sn - 1])
118 w = 0;
119 else {
120 idx = locate_irr(sz, sn, gps->z[ids][iz2]);
121 w =
122 LIN(sz[idx], se[idx], sz[idx + 1], se[idx + 1],
123 gps->z[ids][iz2]);
124 }
125 if (gsl_finite(gps->t[ids][iz2]) && gps->pt[ids][iz2]) {
126 tw += w * gps->t[ids][iz2];
127 wsum += w;
128 gps->pt[ids][iz2] *= w;
129 wmax = GSL_MAX(w, wmax);
130 }
131 }
132 tw /= wsum;
133 for (iz2 = 0; iz2 < gps->nz[ids]; iz2++)
134 gps->pt[ids][iz2] /= wmax;
135 }
136
137 /* Get grid indices... */
138 ix = (int) ((gps->lon[ids][gps->nz[ids] / 2] - lon0)
139 / (lon1 - lon0) * (double) nx);
140 iy = (int) ((gps->lat[ids][gps->nz[ids] / 2] - lat0)
141 / (lat1 - lat0) * (double) ny);
142 if (ix < 0 || ix >= nx || iy < 0 || iy >= ny)
143 continue;
144
145 /* Get mean time and tropopause height... */
146 mtime[ix][iy] += gps->time[ids];
147 thmean[ix][iy] += gps->th[ids];
148 twmean[ix][iy] += tw;
149 np[ix][iy]++;
150
151 /* Loop over altitudes... */
152 for (iz2 = 0; iz2 < gps->nz[ids]; iz2++) {
153
154 /* Get grid indices... */
155 iz = (int) ((gps->z[ids][iz2] - z0)
156 / (z1 - z0) * (double) nz);
157 if (iz < 0 || iz >= nz)
158 continue;
159
160 /* Check data... */
161 if (!gsl_finite(gps->t[ids][iz2]) || !gsl_finite(gps->pt[ids][iz2]))
162 continue;
163
164 /* Get statistics of perturbations... */
165 tmean[ix][iy][iz] += gps->t[ids][iz2];
166 mean[ix][iy][iz] += gps->pt[ids][iz2];
167 var[ix][iy][iz] += gsl_pow_2(gps->pt[ids][iz2]);
168 max[ix][iy][iz] = GSL_MAX(max[ix][iy][iz], gps->pt[ids][iz2]);
169 min[ix][iy][iz] = GSL_MIN(min[ix][iy][iz], gps->pt[ids][iz2]);
170 n[ix][iy][iz]++;
171 }
172 }
173 }
174
175 /* Analyze results... */
176 for (ix = 0; ix < nx; ix++)
177 for (iy = 0; iy < ny; iy++) {
178
179 /* Get mean time and tropopause height... */
180 if (np[ix][iy] > 0) {
181 mtime[ix][iy] /= (double) np[ix][iy];
182 thmean[ix][iy] /= (double) np[ix][iy];
183 twmean[ix][iy] /= (double) np[ix][iy];
184 } else {
185 mtime[ix][iy] = GSL_NAN;
186 thmean[ix][iy] = GSL_NAN;
187 twmean[ix][iy] = GSL_NAN;
188 }
189
190 /* Loop over altitudes... */
191 for (iz = 0; iz < nz; iz++) {
192
193 /* Get geolocation... */
194 gz[iz] = z0 + (iz + 0.5) / (double) nz *(
195 z1 - z0);
196 glon[ix] = lon0 + (ix + 0.5) / (double) nx *(
197 lon1 - lon0);
198 glat[iy] = lat0 + (iy + 0.5) / (double) ny *(
199 lat1 - lat0);
200
201 /* Get mean perturbation and variance... */
202 if (n[ix][iy][iz] > 0) {
203 tmean[ix][iy][iz]
204 /= (double) n[ix][iy][iz];
205 mean[ix][iy][iz]
206 /= (double) n[ix][iy][iz];
207 var[ix][iy][iz]
208 = var[ix][iy][iz] / (double) n[ix][iy][iz]
209 - gsl_pow_2(mean[ix][iy][iz]);
210 } else {
211 tmean[ix][iy][iz] = GSL_NAN;
212 mean[ix][iy][iz] = GSL_NAN;
213 var[ix][iy][iz] = GSL_NAN;
214 min[ix][iy][iz] = GSL_NAN;
215 max[ix][iy][iz] = GSL_NAN;
216 }
217 }
218 }
219
220 /* Create file... */
221 printf("Write variance statistics: %s\n", argv[2]);
222 if (!(out = fopen(argv[2], "w")))
223 ERRMSG("Cannot create file!");
224
225 /* Write header... */
226 fprintf(out,
227 "# $1 = time [s]\n"
228 "# $2 = altitude [km]\n"
229 "# $3 = longitude [deg]\n"
230 "# $4 = latitude [deg]\n"
231 "# $5 = number of profiles\n"
232 "# $6 = number of data points\n"
233 "# $7 = mean perturbation [K]\n"
234 "# $8 = minimum perturbation [K]\n"
235 "# $9 = maximum perturbation [K]\n"
236 "# $10 = variance [K^2]\n"
237 "# $11 = mean temperature [K]\n"
238 "# $12 = mean weighted temperature [K]\n"
239 "# $13 = mean tropopause height [km]\n");
240
241 /* Write results... */
242 for (iz = 0; iz < nz; iz++) {
243 for (iy = 0; iy < ny; iy++) {
244 if (iy == 0 || nx > 1)
245 fprintf(out, "\n");
246 for (ix = 0; ix < nx; ix++)
247 fprintf(out, "%.2f %g %g %g %d %d %g %g %g %g %g %g %g\n",
248 mtime[ix][iy], gz[iz], glon[ix], glat[iy],
249 np[ix][iy], n[ix][iy][iz], mean[ix][iy][iz],
250 min[ix][iy][iz], max[ix][iy][iz], var[ix][iy][iz],
251 tmean[ix][iy][iz], twmean[ix][iy], thmean[ix][iy]);
252 }
253 }
254
255 /* Close file... */
256 fclose(out);
257
258 /* Free... */
259 free(gps);
260
261 return EXIT_SUCCESS;
262}
read_gps
void read_gps(char *filename, gps_t *gps)
Read GPS-RO data file.
Definition: libgps.c:611
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 Sat Jul 19 2025 19:02:12 for GPS Code Collection by doxygen 1.9.4