37 double *ahtd, *aqtd, *avtd, ahtdm, aqtdm[
NQ], avtdm, *lat1_old, *lat2_old,
38 *lh1, *lh2, *lon1_old, *lon2_old, *lv1, *lv2, *rhtd, *rqtd, *rvtd, rhtdm,
39 rqtdm[
NQ], rvtdm, t0 =
40 0, x0[3], x1[3], x2[3], z1, *z1_old, z2, *z2_old, *work;
47 ALLOC(lon1_old,
double,
49 ALLOC(lat1_old,
double,
57 ALLOC(lon2_old,
double,
59 ALLOC(lat2_old,
double,
84 ERRMSG(
"Give parameters: <ctl> <dist.tab> <param> <atm1a> <atm1b>"
85 " [<atm2a> <atm2b> ...]");
90 (int)
scan_ctl(argv[1], argc, argv,
"DIST_ENS", -1,
"-999", NULL);
92 P(
scan_ctl(argv[1], argc, argv,
"DIST_Z0", -1,
"-1000", NULL));
94 P(
scan_ctl(argv[1], argc, argv,
"DIST_Z1", -1,
"1000", NULL));
96 scan_ctl(argv[1], argc, argv,
"DIST_LAT0", -1,
"-1000", NULL);
98 scan_ctl(argv[1], argc, argv,
"DIST_LAT1", -1,
"1000", NULL);
100 scan_ctl(argv[1], argc, argv,
"DIST_LON0", -1,
"-1000", NULL);
102 scan_ctl(argv[1], argc, argv,
"DIST_LON1", -1,
"1000", NULL);
103 const double zscore =
104 scan_ctl(argv[1], argc, argv,
"DIST_ZSCORE", -1,
"-999", NULL);
107 LOG(1,
"Write transport deviations: %s", argv[2]);
110 if (!(out = fopen(argv[2],
"w")))
111 ERRMSG(
"Cannot create file!");
116 "# $2 = time difference [s]\n"
117 "# $3 = absolute horizontal distance (%s) [km]\n"
118 "# $4 = relative horizontal distance (%s) [%%]\n"
119 "# $5 = absolute vertical distance (%s) [km]\n"
120 "# $6 = relative vertical distance (%s) [%%]\n",
121 argv[3], argv[3], argv[3], argv[3]);
122 for (
int iq = 0; iq < ctl.
nq; iq++)
124 "# $%d = %s absolute difference (%s) [%s]\n"
125 "# $%d = %s relative difference (%s) [%%]\n",
127 8 + 2 * iq, ctl.
qnt_name[iq], argv[3]);
128 fprintf(out,
"# $%d = number of particles\n\n", 7 + 2 * ctl.
nq);
131 for (
int f = 4; f < argc; f += 2) {
139 if (atm1->
np != atm2->
np)
140 ERRMSG(
"Different numbers of particles!");
153 for (
int ip = 0; ip < atm1->
np; ip++) {
154 ahtd[ip] = avtd[ip] = rhtd[ip] = rvtd[ip] = 0;
155 for (
int iq = 0; iq < ctl.
nq; iq++)
156 aqtd[iq *
NP + ip] = rqtd[iq *
NP + ip] = 0;
160 for (
int ip = 0; ip < atm1->
np; ip++) {
165 ERRMSG(
"Air parcel index does not match!");
170 || atm2->
q[ctl.
qnt_ens][ip] != ens))
174 if (!isfinite(atm1->
time[ip]) || !isfinite(atm2->
time[ip]))
178 if (atm1->
p[ip] > p0 || atm1->
p[ip] < p1
179 || atm1->
lon[ip] < lon0 || atm1->
lon[ip] > lon1
180 || atm1->
lat[ip] < lat0 || atm1->
lat[ip] > lat1)
182 if (atm2->
p[ip] > p0 || atm2->
p[ip] < p1
183 || atm2->
lon[ip] < lon0 || atm2->
lon[ip] > lon1
184 || atm2->
lat[ip] < lat0 || atm2->
lat[ip] > lat1)
194 ahtd[np] =
DIST(x1, x2);
196 for (
int iq = 0; iq < ctl.
nq; iq++)
197 aqtd[iq *
NP + np] = atm1->
q[iq][ip] - atm2->
q[iq][ip];
203 geo2cart(0, lon1_old[ip], lat1_old[ip], x0);
204 lh1[ip] +=
DIST(x0, x1);
205 lv1[ip] += fabs(z1_old[ip] - z1);
207 geo2cart(0, lon2_old[ip], lat2_old[ip], x0);
208 lh2[ip] +=
DIST(x0, x2);
209 lv2[ip] += fabs(z2_old[ip] - z2);
212 if (lh1[ip] + lh2[ip] > 0)
213 rhtd[np] = 200. *
DIST(x1, x2) / (lh1[ip] + lh2[ip]);
214 if (lv1[ip] + lv2[ip] > 0)
215 rvtd[np] = 200. * (z1 - z2) / (lv1[ip] + lv2[ip]);
219 for (
int iq = 0; iq < ctl.
nq; iq++)
220 rqtd[iq *
NP + np] = 200. * (atm1->
q[iq][ip] - atm2->
q[iq][ip])
221 / (fabs(atm1->
q[iq][ip]) + fabs(atm2->
q[iq][ip]));
224 lon1_old[ip] = atm1->
lon[ip];
225 lat1_old[ip] = atm1->
lat[ip];
228 lon2_old[ip] = atm2->
lon[ip];
229 lat2_old[ip] = atm2->
lat[ip];
237 if (zscore > 0 && np > 1) {
240 const size_t n = (size_t) np;
241 const double muh = gsl_stats_mean(ahtd, 1, n);
242 const double muv = gsl_stats_mean(avtd, 1, n);
243 const double sigh = gsl_stats_sd(ahtd, 1, n);
244 const double sigv = gsl_stats_sd(avtd, 1, n);
248 for (
size_t i = 0; i < n; i++)
249 if (fabs((ahtd[i] - muh) / sigh) < zscore
250 && fabs((avtd[i] - muv) / sigv) < zscore) {
255 for (
int iq = 0; iq < ctl.
nq; iq++) {
256 aqtd[iq *
NP + np] = aqtd[iq *
NP + (int) i];
257 rqtd[iq *
NP + np] = rqtd[iq *
NP + (int) i];
264 if (strcasecmp(argv[3],
"mean") == 0) {
265 ahtdm = gsl_stats_mean(ahtd, 1, (
size_t) np);
266 rhtdm = gsl_stats_mean(rhtd, 1, (
size_t) np);
267 avtdm = gsl_stats_mean(avtd, 1, (
size_t) np);
268 rvtdm = gsl_stats_mean(rvtd, 1, (
size_t) np);
269 for (
int iq = 0; iq < ctl.
nq; iq++) {
270 aqtdm[iq] = gsl_stats_mean(&aqtd[iq *
NP], 1, (
size_t) np);
271 rqtdm[iq] = gsl_stats_mean(&rqtd[iq *
NP], 1, (
size_t) np);
273 }
else if (strcasecmp(argv[3],
"stddev") == 0) {
274 ahtdm = gsl_stats_sd(ahtd, 1, (
size_t) np);
275 rhtdm = gsl_stats_sd(rhtd, 1, (
size_t) np);
276 avtdm = gsl_stats_sd(avtd, 1, (
size_t) np);
277 rvtdm = gsl_stats_sd(rvtd, 1, (
size_t) np);
278 for (
int iq = 0; iq < ctl.
nq; iq++) {
279 aqtdm[iq] = gsl_stats_sd(&aqtd[iq *
NP], 1, (
size_t) np);
280 rqtdm[iq] = gsl_stats_sd(&rqtd[iq *
NP], 1, (
size_t) np);
282 }
else if (strcasecmp(argv[3],
"min") == 0) {
283 ahtdm = gsl_stats_min(ahtd, 1, (
size_t) np);
284 rhtdm = gsl_stats_min(rhtd, 1, (
size_t) np);
285 avtdm = gsl_stats_min(avtd, 1, (
size_t) np);
286 rvtdm = gsl_stats_min(rvtd, 1, (
size_t) np);
287 for (
int iq = 0; iq < ctl.
nq; iq++) {
288 aqtdm[iq] = gsl_stats_min(&aqtd[iq *
NP], 1, (
size_t) np);
289 rqtdm[iq] = gsl_stats_min(&rqtd[iq *
NP], 1, (
size_t) np);
291 }
else if (strcasecmp(argv[3],
"max") == 0) {
292 ahtdm = gsl_stats_max(ahtd, 1, (
size_t) np);
293 rhtdm = gsl_stats_max(rhtd, 1, (
size_t) np);
294 avtdm = gsl_stats_max(avtd, 1, (
size_t) np);
295 rvtdm = gsl_stats_max(rvtd, 1, (
size_t) np);
296 for (
int iq = 0; iq < ctl.
nq; iq++) {
297 aqtdm[iq] = gsl_stats_max(&aqtd[iq *
NP], 1, (
size_t) np);
298 rqtdm[iq] = gsl_stats_max(&rqtd[iq *
NP], 1, (
size_t) np);
300 }
else if (strcasecmp(argv[3],
"skew") == 0) {
301 ahtdm = gsl_stats_skew(ahtd, 1, (
size_t) np);
302 rhtdm = gsl_stats_skew(rhtd, 1, (
size_t) np);
303 avtdm = gsl_stats_skew(avtd, 1, (
size_t) np);
304 rvtdm = gsl_stats_skew(rvtd, 1, (
size_t) np);
305 for (
int iq = 0; iq < ctl.
nq; iq++) {
306 aqtdm[iq] = gsl_stats_skew(&aqtd[iq *
NP], 1, (
size_t) np);
307 rqtdm[iq] = gsl_stats_skew(&rqtd[iq *
NP], 1, (
size_t) np);
309 }
else if (strcasecmp(argv[3],
"kurt") == 0) {
310 ahtdm = gsl_stats_kurtosis(ahtd, 1, (
size_t) np);
311 rhtdm = gsl_stats_kurtosis(rhtd, 1, (
size_t) np);
312 avtdm = gsl_stats_kurtosis(avtd, 1, (
size_t) np);
313 rvtdm = gsl_stats_kurtosis(rvtd, 1, (
size_t) np);
314 for (
int iq = 0; iq < ctl.
nq; iq++) {
315 aqtdm[iq] = gsl_stats_kurtosis(&aqtd[iq *
NP], 1, (
size_t) np);
316 rqtdm[iq] = gsl_stats_kurtosis(&rqtd[iq *
NP], 1, (
size_t) np);
318 }
else if (strcasecmp(argv[3],
"absdev") == 0) {
319 ahtdm = gsl_stats_absdev_m(ahtd, 1, (
size_t) np, 0.0);
320 rhtdm = gsl_stats_absdev_m(rhtd, 1, (
size_t) np, 0.0);
321 avtdm = gsl_stats_absdev_m(avtd, 1, (
size_t) np, 0.0);
322 rvtdm = gsl_stats_absdev_m(rvtd, 1, (
size_t) np, 0.0);
323 for (
int iq = 0; iq < ctl.
nq; iq++) {
324 aqtdm[iq] = gsl_stats_absdev_m(&aqtd[iq *
NP], 1, (
size_t) np, 0.0);
325 rqtdm[iq] = gsl_stats_absdev_m(&rqtd[iq *
NP], 1, (
size_t) np, 0.0);
327 }
else if (strcasecmp(argv[3],
"median") == 0) {
328 ahtdm = gsl_stats_median(ahtd, 1, (
size_t) np);
329 rhtdm = gsl_stats_median(rhtd, 1, (
size_t) np);
330 avtdm = gsl_stats_median(avtd, 1, (
size_t) np);
331 rvtdm = gsl_stats_median(rvtd, 1, (
size_t) np);
332 for (
int iq = 0; iq < ctl.
nq; iq++) {
333 aqtdm[iq] = gsl_stats_median(&aqtd[iq *
NP], 1, (
size_t) np);
334 rqtdm[iq] = gsl_stats_median(&rqtd[iq *
NP], 1, (
size_t) np);
336 }
else if (strcasecmp(argv[3],
"mad") == 0) {
337 ahtdm = gsl_stats_mad0(ahtd, 1, (
size_t) np, work);
338 rhtdm = gsl_stats_mad0(rhtd, 1, (
size_t) np, work);
339 avtdm = gsl_stats_mad0(avtd, 1, (
size_t) np, work);
340 rvtdm = gsl_stats_mad0(rvtd, 1, (
size_t) np, work);
341 for (
int iq = 0; iq < ctl.
nq; iq++) {
342 aqtdm[iq] = gsl_stats_mad0(&aqtd[iq *
NP], 1, (
size_t) np, work);
343 rqtdm[iq] = gsl_stats_mad0(&rqtd[iq *
NP], 1, (
size_t) np, work);
346 ERRMSG(
"Unknown parameter!");
349 fprintf(out,
"%.2f %.2f %g %g %g %g", t, t - t0,
350 ahtdm, rhtdm, avtdm, rvtdm);
351 for (
int iq = 0; iq < ctl.
nq; iq++) {
357 fprintf(out,
" %d\n", np);
int main(int argc, char *argv[])
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.
MPTRAC library declarations.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
#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.