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;
68 ERRMSG(
"Missing or invalid command-line arguments.\n\n"
69 "Usage: atm_dist <ctl> <dist.tab> <param> <atm1a> <atm1b> [<atm2a> <atm2b> ...]\n\n"
70 "Use -h for full help.");
75 ALLOC(lon1_old,
double,
77 ALLOC(lat1_old,
double,
85 ALLOC(lon2_old,
double,
87 ALLOC(lat2_old,
double,
114 ERRMSG(
"atm_dist currently supports only lat/lon grids");
116 (int)
scan_ctl(argv[1], argc, argv,
"DIST_ENS", -1,
"-999", NULL);
118 P(
scan_ctl(argv[1], argc, argv,
"DIST_Z0", -1,
"-1000", NULL));
120 P(
scan_ctl(argv[1], argc, argv,
"DIST_Z1", -1,
"1000", NULL));
122 scan_ctl(argv[1], argc, argv,
"DIST_LAT0", -1,
"-1000", NULL);
124 scan_ctl(argv[1], argc, argv,
"DIST_LAT1", -1,
"1000", NULL);
126 scan_ctl(argv[1], argc, argv,
"DIST_LON0", -1,
"-1000", NULL);
128 scan_ctl(argv[1], argc, argv,
"DIST_LON1", -1,
"1000", NULL);
129 const double zscore =
130 scan_ctl(argv[1], argc, argv,
"DIST_ZSCORE", -1,
"-999", NULL);
133 LOG(1,
"Write transport deviations: %s", argv[2]);
136 if (!(out = fopen(argv[2],
"w")))
137 ERRMSG(
"Cannot create file!");
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++) {
150 scan_ctl(argv[1], argc, argv,
"DIST_REL_MIN", iq,
"0", NULL);
153 "# Relative %s differences are masked where |q1| + |q2| <= %g %s.\n",
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]);
161 fprintf(out,
"# $%d = number of particles\n\n", 7 + 2 * ctl.
nq);
164 for (
int f = 4; f < argc; f += 2) {
172 if (atm1->
np != atm2->
np)
173 ERRMSG(
"Different numbers of particles!");
176 int time_offset = ctl.
atm_type < 2 ? 23 : 22;
187 for (
int ip = 0; ip < atm1->
np; ip++) {
188 ahtd[ip] = avtd[ip] = rhtd[ip] = rvtd[ip] = 0;
189 for (
int iq = 0; iq < ctl.
nq; iq++)
190 aqtd[iq *
NP + ip] = rqtd[iq *
NP + ip] = 0;
194 for (
int ip = 0; ip < atm1->
np; ip++) {
199 ERRMSG(
"Air parcel index does not match!");
204 || atm2->
q[ctl.
qnt_ens][ip] != ens))
208 if (!isfinite(atm1->
time[ip]) || !isfinite(atm2->
time[ip]))
212 if (atm1->
p[ip] > p0 || atm1->
p[ip] < p1
213 || atm1->
lon[ip] < lon0 || atm1->
lon[ip] > lon1
214 || atm1->
lat[ip] < lat0 || atm1->
lat[ip] > lat1)
216 if (atm2->
p[ip] > p0 || atm2->
p[ip] < p1
217 || atm2->
lon[ip] < lon0 || atm2->
lon[ip] > lon1
218 || atm2->
lat[ip] < lat0 || atm2->
lat[ip] > lat1)
228 ahtd[np] =
DIST(x1, x2);
230 for (
int iq = 0; iq < ctl.
nq; iq++)
231 aqtd[iq *
NP + np] = atm1->
q[iq][ip] - atm2->
q[iq][ip];
237 geo2cart(0, lon1_old[ip], lat1_old[ip], x0);
238 lh1[ip] +=
DIST(x0, x1);
239 lv1[ip] += fabs(z1_old[ip] - z1);
241 geo2cart(0, lon2_old[ip], lat2_old[ip], x0);
242 lh2[ip] +=
DIST(x0, x2);
243 lv2[ip] += fabs(z2_old[ip] - z2);
246 if (lh1[ip] + lh2[ip] > 0)
247 rhtd[np] = 200. *
DIST(x1, x2) / (lh1[ip] + lh2[ip]);
248 if (lv1[ip] + lv2[ip] > 0)
249 rvtd[np] = 200. * (z1 - z2) / (lv1[ip] + lv2[ip]);
253 for (
int iq = 0; iq < ctl.
nq; iq++) {
254 const double q1 = atm1->
q[iq][ip];
255 const double q2 = atm2->
q[iq][ip];
256 const double denom = fabs(q1) + fabs(q2);
257 if (denom <= rel_min[iq])
258 rqtd[iq *
NP + np] = NAN;
260 rqtd[iq *
NP + np] = 200. * (q1 - q2) / denom;
264 lon1_old[ip] = atm1->
lon[ip];
265 lat1_old[ip] = atm1->
lat[ip];
268 lon2_old[ip] = atm2->
lon[ip];
269 lat2_old[ip] = atm2->
lat[ip];
277 if (zscore > 0 && np > 1) {
280 const size_t n = (size_t) np;
281 const double muh = gsl_stats_mean(ahtd, 1, n);
282 const double muv = gsl_stats_mean(avtd, 1, n);
283 const double sigh = gsl_stats_sd(ahtd, 1, n);
284 const double sigv = gsl_stats_sd(avtd, 1, n);
288 for (
size_t i = 0; i < n; i++)
289 if (fabs((ahtd[i] - muh) / sigh) < zscore
290 && fabs((avtd[i] - muv) / sigv) < zscore) {
295 for (
int iq = 0; iq < ctl.
nq; iq++) {
296 aqtd[iq *
NP + np] = aqtd[iq *
NP + (int) i];
297 rqtd[iq *
NP + np] = rqtd[iq *
NP + (int) i];
304 if (strcasecmp(argv[3],
"mean") == 0) {
305 ahtdm = gsl_stats_mean(ahtd, 1, (
size_t) np);
306 rhtdm = gsl_stats_mean(rhtd, 1, (
size_t) np);
307 avtdm = gsl_stats_mean(avtd, 1, (
size_t) np);
308 rvtdm = gsl_stats_mean(rvtd, 1, (
size_t) np);
309 for (
int iq = 0; iq < ctl.
nq; iq++) {
310 aqtdm[iq] = gsl_stats_mean(&aqtd[iq *
NP], 1, (
size_t) np);
313 }
else if (strcasecmp(argv[3],
"stddev") == 0) {
314 ahtdm = gsl_stats_sd(ahtd, 1, (
size_t) np);
315 rhtdm = gsl_stats_sd(rhtd, 1, (
size_t) np);
316 avtdm = gsl_stats_sd(avtd, 1, (
size_t) np);
317 rvtdm = gsl_stats_sd(rvtd, 1, (
size_t) np);
318 for (
int iq = 0; iq < ctl.
nq; iq++) {
319 aqtdm[iq] = gsl_stats_sd(&aqtd[iq *
NP], 1, (
size_t) np);
322 }
else if (strcasecmp(argv[3],
"min") == 0) {
323 ahtdm = gsl_stats_min(ahtd, 1, (
size_t) np);
324 rhtdm = gsl_stats_min(rhtd, 1, (
size_t) np);
325 avtdm = gsl_stats_min(avtd, 1, (
size_t) np);
326 rvtdm = gsl_stats_min(rvtd, 1, (
size_t) np);
327 for (
int iq = 0; iq < ctl.
nq; iq++) {
328 aqtdm[iq] = gsl_stats_min(&aqtd[iq *
NP], 1, (
size_t) np);
331 }
else if (strcasecmp(argv[3],
"max") == 0) {
332 ahtdm = gsl_stats_max(ahtd, 1, (
size_t) np);
333 rhtdm = gsl_stats_max(rhtd, 1, (
size_t) np);
334 avtdm = gsl_stats_max(avtd, 1, (
size_t) np);
335 rvtdm = gsl_stats_max(rvtd, 1, (
size_t) np);
336 for (
int iq = 0; iq < ctl.
nq; iq++) {
337 aqtdm[iq] = gsl_stats_max(&aqtd[iq *
NP], 1, (
size_t) np);
340 }
else if (strcasecmp(argv[3],
"skew") == 0) {
341 ahtdm = gsl_stats_skew(ahtd, 1, (
size_t) np);
342 rhtdm = gsl_stats_skew(rhtd, 1, (
size_t) np);
343 avtdm = gsl_stats_skew(avtd, 1, (
size_t) np);
344 rvtdm = gsl_stats_skew(rvtd, 1, (
size_t) np);
345 for (
int iq = 0; iq < ctl.
nq; iq++) {
346 aqtdm[iq] = gsl_stats_skew(&aqtd[iq *
NP], 1, (
size_t) np);
349 }
else if (strcasecmp(argv[3],
"kurt") == 0) {
350 ahtdm = gsl_stats_kurtosis(ahtd, 1, (
size_t) np);
351 rhtdm = gsl_stats_kurtosis(rhtd, 1, (
size_t) np);
352 avtdm = gsl_stats_kurtosis(avtd, 1, (
size_t) np);
353 rvtdm = gsl_stats_kurtosis(rvtd, 1, (
size_t) np);
354 for (
int iq = 0; iq < ctl.
nq; iq++) {
355 aqtdm[iq] = gsl_stats_kurtosis(&aqtd[iq *
NP], 1, (
size_t) np);
358 }
else if (strcasecmp(argv[3],
"absdev") == 0) {
359 ahtdm = gsl_stats_absdev_m(ahtd, 1, (
size_t) np, 0.0);
360 rhtdm = gsl_stats_absdev_m(rhtd, 1, (
size_t) np, 0.0);
361 avtdm = gsl_stats_absdev_m(avtd, 1, (
size_t) np, 0.0);
362 rvtdm = gsl_stats_absdev_m(rvtd, 1, (
size_t) np, 0.0);
363 for (
int iq = 0; iq < ctl.
nq; iq++) {
364 aqtdm[iq] = gsl_stats_absdev_m(&aqtd[iq *
NP], 1, (
size_t) np, 0.0);
367 }
else if (strcasecmp(argv[3],
"median") == 0) {
368 ahtdm = gsl_stats_median(ahtd, 1, (
size_t) np);
369 rhtdm = gsl_stats_median(rhtd, 1, (
size_t) np);
370 avtdm = gsl_stats_median(avtd, 1, (
size_t) np);
371 rvtdm = gsl_stats_median(rvtd, 1, (
size_t) np);
372 for (
int iq = 0; iq < ctl.
nq; iq++) {
373 aqtdm[iq] = gsl_stats_median(&aqtd[iq *
NP], 1, (
size_t) np);
376 }
else if (strcasecmp(argv[3],
"mad") == 0) {
377 ahtdm = gsl_stats_mad0(ahtd, 1, (
size_t) np, work);
378 rhtdm = gsl_stats_mad0(rhtd, 1, (
size_t) np, work);
379 avtdm = gsl_stats_mad0(avtd, 1, (
size_t) np, work);
380 rvtdm = gsl_stats_mad0(rvtd, 1, (
size_t) np, work);
381 for (
int iq = 0; iq < ctl.
nq; iq++) {
382 aqtdm[iq] = gsl_stats_mad0(&aqtd[iq *
NP], 1, (
size_t) np, work);
386 ERRMSG(
"Unknown parameter!");
389 fprintf(out,
"%.2f %.2f %g %g %g %g", t, t - t0,
390 ahtdm, rhtdm, avtdm, rvtdm);
391 for (
int iq = 0; iq < ctl.
nq; iq++) {
397 fprintf(out,
" %d\n", np);