48 {
49
51
53
54 FILE *out;
55
56 double *ahtd, *aqtd, *avtd, ahtdm, aqtdm[
NQ], avtdm, *lat1_old, *lat2_old,
57 *lh1, *lh2, *lon1_old, *lon2_old, *lv1, *lv2, *rhtd, *rqtd, *rvtd, rhtdm,
58 rel_min[
NQ], rqtdm[
NQ], rvtdm, t0 =
59 0, x0[3], x1[3], x2[3], z1, *z1_old, z2, *z2_old, *work;
60
61 int init = 0, np;
62
63
66 ALLOC(lon1_old,
double,
68 ALLOC(lat1_old,
double,
76 ALLOC(lon2_old,
double,
78 ALLOC(lat2_old,
double,
100
101
103
104
105 if (argc < 6)
106 ERRMSG(
"Missing or invalid command-line arguments.\n\n"
107 "Usage: atm_dist <ctl> <dist.tab> <param> <atm1a> <atm1b> [<atm2a> <atm2b> ...]\n\n"
108 "Use -h for full help.");
109
110
112 const int ens =
113 (int)
scan_ctl(argv[1], argc, argv,
"DIST_ENS", -1,
"-999", NULL);
114 const double p0 =
115 P(
scan_ctl(argv[1], argc, argv,
"DIST_Z0", -1,
"-1000", NULL));
116 const double p1 =
117 P(
scan_ctl(argv[1], argc, argv,
"DIST_Z1", -1,
"1000", NULL));
118 const double lat0 =
119 scan_ctl(argv[1], argc, argv,
"DIST_LAT0", -1,
"-1000", NULL);
120 const double lat1 =
121 scan_ctl(argv[1], argc, argv,
"DIST_LAT1", -1,
"1000", NULL);
122 const double lon0 =
123 scan_ctl(argv[1], argc, argv,
"DIST_LON0", -1,
"-1000", NULL);
124 const double lon1 =
125 scan_ctl(argv[1], argc, argv,
"DIST_LON1", -1,
"1000", NULL);
126 const double zscore =
127 scan_ctl(argv[1], argc, argv,
"DIST_ZSCORE", -1,
"-999", NULL);
128
129 for (
int iq = 0; iq <
NQ; iq++)
130 rel_min[iq] = 0;
131
132
133 LOG(1,
"Write transport deviations: %s", argv[2]);
134
135
136 if (!(out = fopen(argv[2], "w")))
137 ERRMSG(
"Cannot create file!");
138
139
140 fprintf(out,
141 "# $1 = time [s]\n"
142 "# $2 = time difference [s]\n"
143 "# $3 = absolute horizontal distance (%s) [km]\n"
144 "# $4 = relative horizontal distance (%s) [%%]\n"
145 "# $5 = absolute vertical distance (%s) [km]\n"
146 "# $6 = relative vertical distance (%s) [%%]\n",
147 argv[3], argv[3], argv[3], argv[3]);
148 for (
int iq = 0; iq < ctl.
nq; iq++) {
149 rel_min[iq] =
150 scan_ctl(argv[1], argc, argv,
"DIST_REL_MIN", iq,
"0", NULL);
151 if (rel_min[iq] > 0)
152 fprintf(out,
153 "# Relative %s differences are masked where |q1| + |q2| <= %g %s.\n",
155 fprintf(out,
156 "# $%d = %s absolute difference (%s) [%s]\n"
157 "# $%d = %s relative difference (%s) [%%]\n",
159 8 + 2 * iq, ctl.
qnt_name[iq], argv[3]);
160 }
161 fprintf(out,
"# $%d = number of particles\n\n", 7 + 2 * ctl.
nq);
162
163
164 for (int f = 4; f < argc; f += 2) {
165
166
169 continue;
170
171
172 if (atm1->
np != atm2->
np)
173 ERRMSG(
"Different numbers of particles!");
174
175
177
178
179 if (!init) {
180 init = 1;
181 t0 = t;
182 }
183
184
185 np = 0;
186 for (
int ip = 0; ip < atm1->
np; ip++) {
187 ahtd[ip] = avtd[ip] = rhtd[ip] = rvtd[ip] = 0;
188 for (
int iq = 0; iq < ctl.
nq; iq++)
189 aqtd[iq *
NP + ip] = rqtd[iq *
NP + ip] = 0;
190 }
191
192
193 for (
int ip = 0; ip < atm1->
np; ip++) {
194
195
198 ERRMSG(
"Air parcel index does not match!");
199
200
203 || atm2->
q[ctl.
qnt_ens][ip] != ens))
204 continue;
205
206
207 if (!isfinite(atm1->
time[ip]) || !isfinite(atm2->
time[ip]))
208 continue;
209
210
211 if (atm1->
p[ip] > p0 || atm1->
p[ip] < p1
212 || atm1->
lon[ip] < lon0 || atm1->
lon[ip] > lon1
213 || atm1->
lat[ip] < lat0 || atm1->
lat[ip] > lat1)
214 continue;
215 if (atm2->
p[ip] > p0 || atm2->
p[ip] < p1
216 || atm2->
lon[ip] < lon0 || atm2->
lon[ip] > lon1
217 || atm2->
lat[ip] < lat0 || atm2->
lat[ip] > lat1)
218 continue;
219
220
225
226
227 ahtd[np] =
DIST(x1, x2);
228 avtd[np] = z1 - z2;
229 for (
int iq = 0; iq < ctl.
nq; iq++)
230 aqtd[iq *
NP + np] = atm1->
q[iq][ip] - atm2->
q[iq][ip];
231
232
233 if (f > 4) {
234
235
236 geo2cart(0, lon1_old[ip], lat1_old[ip], x0);
237 lh1[ip] +=
DIST(x0, x1);
238 lv1[ip] += fabs(z1_old[ip] - z1);
239
240 geo2cart(0, lon2_old[ip], lat2_old[ip], x0);
241 lh2[ip] +=
DIST(x0, x2);
242 lv2[ip] += fabs(z2_old[ip] - z2);
243
244
245 if (lh1[ip] + lh2[ip] > 0)
246 rhtd[np] = 200. *
DIST(x1, x2) / (lh1[ip] + lh2[ip]);
247 if (lv1[ip] + lv2[ip] > 0)
248 rvtd[np] = 200. * (z1 - z2) / (lv1[ip] + lv2[ip]);
249 }
250
251
252 for (
int iq = 0; iq < ctl.
nq; iq++) {
253 const double q1 = atm1->
q[iq][ip];
254 const double q2 = atm2->
q[iq][ip];
255 const double denom = fabs(q1) + fabs(q2);
256 if (denom <= rel_min[iq])
257 rqtd[iq *
NP + np] = GSL_NAN;
258 else
259 rqtd[iq *
NP + np] = 200. * (q1 - q2) / denom;
260 }
261
262
263 lon1_old[ip] = atm1->
lon[ip];
264 lat1_old[ip] = atm1->
lat[ip];
265 z1_old[ip] = z1;
266
267 lon2_old[ip] = atm2->
lon[ip];
268 lat2_old[ip] = atm2->
lat[ip];
269 z2_old[ip] = z2;
270
271
272 np++;
273 }
274
275
276 if (zscore > 0 && np > 1) {
277
278
279 const size_t n = (size_t) np;
280 const double muh = gsl_stats_mean(ahtd, 1, n);
281 const double muv = gsl_stats_mean(avtd, 1, n);
282 const double sigh = gsl_stats_sd(ahtd, 1, n);
283 const double sigv = gsl_stats_sd(avtd, 1, n);
284
285
286 np = 0;
287 for (size_t i = 0; i < n; i++)
288 if (fabs((ahtd[i] - muh) / sigh) < zscore
289 && fabs((avtd[i] - muv) / sigv) < zscore) {
290 ahtd[np] = ahtd[i];
291 rhtd[np] = rhtd[i];
292 avtd[np] = avtd[i];
293 rvtd[np] = rvtd[i];
294 for (
int iq = 0; iq < ctl.
nq; iq++) {
295 aqtd[iq *
NP + np] = aqtd[iq *
NP + (int) i];
296 rqtd[iq *
NP + np] = rqtd[iq *
NP + (int) i];
297 }
298 np++;
299 }
300 }
301
302
303 if (strcasecmp(argv[3], "mean") == 0) {
304 ahtdm = gsl_stats_mean(ahtd, 1, (size_t) np);
305 rhtdm = gsl_stats_mean(rhtd, 1, (size_t) np);
306 avtdm = gsl_stats_mean(avtd, 1, (size_t) np);
307 rvtdm = gsl_stats_mean(rvtd, 1, (size_t) np);
308 for (
int iq = 0; iq < ctl.
nq; iq++) {
309 aqtdm[iq] = gsl_stats_mean(&aqtd[iq *
NP], 1, (
size_t) np);
311 }
312 } else if (strcasecmp(argv[3], "stddev") == 0) {
313 ahtdm = gsl_stats_sd(ahtd, 1, (size_t) np);
314 rhtdm = gsl_stats_sd(rhtd, 1, (size_t) np);
315 avtdm = gsl_stats_sd(avtd, 1, (size_t) np);
316 rvtdm = gsl_stats_sd(rvtd, 1, (size_t) np);
317 for (
int iq = 0; iq < ctl.
nq; iq++) {
318 aqtdm[iq] = gsl_stats_sd(&aqtd[iq *
NP], 1, (
size_t) np);
320 }
321 } else if (strcasecmp(argv[3], "min") == 0) {
322 ahtdm = gsl_stats_min(ahtd, 1, (size_t) np);
323 rhtdm = gsl_stats_min(rhtd, 1, (size_t) np);
324 avtdm = gsl_stats_min(avtd, 1, (size_t) np);
325 rvtdm = gsl_stats_min(rvtd, 1, (size_t) np);
326 for (
int iq = 0; iq < ctl.
nq; iq++) {
327 aqtdm[iq] = gsl_stats_min(&aqtd[iq *
NP], 1, (
size_t) np);
329 }
330 } else if (strcasecmp(argv[3], "max") == 0) {
331 ahtdm = gsl_stats_max(ahtd, 1, (size_t) np);
332 rhtdm = gsl_stats_max(rhtd, 1, (size_t) np);
333 avtdm = gsl_stats_max(avtd, 1, (size_t) np);
334 rvtdm = gsl_stats_max(rvtd, 1, (size_t) np);
335 for (
int iq = 0; iq < ctl.
nq; iq++) {
336 aqtdm[iq] = gsl_stats_max(&aqtd[iq *
NP], 1, (
size_t) np);
338 }
339 } else if (strcasecmp(argv[3], "skew") == 0) {
340 ahtdm = gsl_stats_skew(ahtd, 1, (size_t) np);
341 rhtdm = gsl_stats_skew(rhtd, 1, (size_t) np);
342 avtdm = gsl_stats_skew(avtd, 1, (size_t) np);
343 rvtdm = gsl_stats_skew(rvtd, 1, (size_t) np);
344 for (
int iq = 0; iq < ctl.
nq; iq++) {
345 aqtdm[iq] = gsl_stats_skew(&aqtd[iq *
NP], 1, (
size_t) np);
347 }
348 } else if (strcasecmp(argv[3], "kurt") == 0) {
349 ahtdm = gsl_stats_kurtosis(ahtd, 1, (size_t) np);
350 rhtdm = gsl_stats_kurtosis(rhtd, 1, (size_t) np);
351 avtdm = gsl_stats_kurtosis(avtd, 1, (size_t) np);
352 rvtdm = gsl_stats_kurtosis(rvtd, 1, (size_t) np);
353 for (
int iq = 0; iq < ctl.
nq; iq++) {
354 aqtdm[iq] = gsl_stats_kurtosis(&aqtd[iq *
NP], 1, (
size_t) np);
356 }
357 } else if (strcasecmp(argv[3], "absdev") == 0) {
358 ahtdm = gsl_stats_absdev_m(ahtd, 1, (size_t) np, 0.0);
359 rhtdm = gsl_stats_absdev_m(rhtd, 1, (size_t) np, 0.0);
360 avtdm = gsl_stats_absdev_m(avtd, 1, (size_t) np, 0.0);
361 rvtdm = gsl_stats_absdev_m(rvtd, 1, (size_t) np, 0.0);
362 for (
int iq = 0; iq < ctl.
nq; iq++) {
363 aqtdm[iq] = gsl_stats_absdev_m(&aqtd[iq *
NP], 1, (
size_t) np, 0.0);
365 }
366 } else if (strcasecmp(argv[3], "median") == 0) {
367 ahtdm = gsl_stats_median(ahtd, 1, (size_t) np);
368 rhtdm = gsl_stats_median(rhtd, 1, (size_t) np);
369 avtdm = gsl_stats_median(avtd, 1, (size_t) np);
370 rvtdm = gsl_stats_median(rvtd, 1, (size_t) np);
371 for (
int iq = 0; iq < ctl.
nq; iq++) {
372 aqtdm[iq] = gsl_stats_median(&aqtd[iq *
NP], 1, (
size_t) np);
374 }
375 } else if (strcasecmp(argv[3], "mad") == 0) {
376 ahtdm = gsl_stats_mad0(ahtd, 1, (size_t) np, work);
377 rhtdm = gsl_stats_mad0(rhtd, 1, (size_t) np, work);
378 avtdm = gsl_stats_mad0(avtd, 1, (size_t) np, work);
379 rvtdm = gsl_stats_mad0(rvtd, 1, (size_t) np, work);
380 for (
int iq = 0; iq < ctl.
nq; iq++) {
381 aqtdm[iq] = gsl_stats_mad0(&aqtd[iq *
NP], 1, (
size_t) np, work);
383 }
384 } else
385 ERRMSG(
"Unknown parameter!");
386
387
388 fprintf(out, "%.2f %.2f %g %g %g %g", t, t - t0,
389 ahtdm, rhtdm, avtdm, rvtdm);
390 for (
int iq = 0; iq < ctl.
nq; iq++) {
391 fprintf(out, " ");
393 fprintf(out, " ");
395 }
396 fprintf(out, " %d\n", np);
397 }
398
399
400 fclose(out);
401
402
403 free(atm1);
404 free(atm2);
405 free(lon1_old);
406 free(lat1_old);
407 free(z1_old);
408 free(lh1);
409 free(lv1);
410 free(lon2_old);
411 free(lat2_old);
412 free(z2_old);
413 free(lh2);
414 free(lv2);
415 free(ahtd);
416 free(avtd);
417 free(aqtd);
418 free(rhtd);
419 free(rvtd);
420 free(rqtd);
421 free(work);
422
423 return EXIT_SUCCESS;
424}
double finite_stat(const double *data, int np, const char *param, double *work)
Calculate a statistic on finite values only.
double scan_ctl(const char *filename, int argc, char *argv[], const char *varname, const int arridx, const char *defvalue, char *value)
Scans a control file or command-line arguments for a specified variable.
double time_from_filename(const char *filename, const int offset)
Extracts and converts a timestamp from a filename to Julian seconds.
int mptrac_read_atm(const char *filename, const ctl_t *ctl, atm_t *atm)
Reads air parcel data from a specified file into the given atmospheric structure.
void mptrac_read_ctl(const char *filename, int argc, char *argv[], ctl_t *ctl)
Reads control parameters from a configuration file and populates the given structure.
void geo2cart(const double z, const double lon, const double lat, double *x)
Converts geographic coordinates (longitude, latitude, altitude) to Cartesian coordinates.
#define USAGE
Print usage information on -h or --help.
#define Z(p)
Convert pressure to altitude.
#define P(z)
Compute pressure at given altitude.
#define NQ
Maximum number of quantities per data point.
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
#define NP
Maximum number of atmospheric data points.
#define LOG(level,...)
Print a log message with a specified logging level.
#define DIST(a, b)
Calculate the distance between two points in Cartesian coordinates.
double lat[NP]
Latitude [deg].
double lon[NP]
Longitude [deg].
int np
Number of air parcels.
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
double p[NP]
Pressure [hPa].
char qnt_format[NQ][LEN]
Quantity output format.
int atm_type
Type of atmospheric data files (0=ASCII, 1=binary, 2=netCDF, 3=CLaMS_traj, 4=CLaMS_pos).
char qnt_unit[NQ][LEN]
Quantity units.
char qnt_name[NQ][LEN]
Quantity names.
int qnt_ens
Quantity array index for ensemble IDs.
int qnt_idx
Quantity array index for air parcel IDs.
int nq
Number of quantities.