MPTRAC
Functions
atm_stat.c File Reference

Calculate air parcel statistics. More...

#include "mptrac.h"

Go to the source code of this file.

Functions

void usage (void)
 Print command-line help. More...
 
int main (int argc, char *argv[])
 

Detailed Description

Calculate air parcel statistics.

Definition in file atm_stat.c.

Function Documentation

◆ usage()

void usage ( void  )

Print command-line help.

Definition at line 239 of file atm_stat.c.

240 {
241
242 printf("\nMPTRAC atm_stat tool.\n\n");
243 printf("Calculate air parcel statistics.\n");
244 printf("\n");
245 printf("Usage:\n");
246 printf(" atm_stat <ctl> <stat.tab> <param> <atm1> [<atm2> ...]\n");
247 printf("\n");
248 printf("Arguments:\n");
249 printf(" <ctl> Control file.\n");
250 printf(" <stat.tab> Output table.\n");
251 printf
252 (" <param> Statistic: mean, stddev, min, max, skew, kurt, absdev,\n");
253 printf(" median, or mad.\n");
254 printf(" <atm*> Atmospheric input files.\n");
255 printf("\nFurther information:\n");
256 printf(" Manual: https://slcs-jsc.github.io/mptrac/\n");
257}

◆ main()

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

Definition at line 39 of file atm_stat.c.

41 {
42
43 ctl_t ctl;
44
45 atm_t *atm, *atm_filt;
46
47 FILE *out;
48
49 double latm, lonm, t, t0 = NAN, qm[NQ], *work, zm, *zs;
50
51 int init = 0;
52
53 /* Print usage information... */
54 USAGE;
55
56 /* Check arguments... */
57 if (argc < 4)
58 ERRMSG("Missing or invalid command-line arguments.\n\n"
59 "Usage: atm_stat <ctl> <stat.tab> <param> <atm1> [<atm2> ...]\n\n"
60 "Use -h for full help.");
61
62 /* Allocate... */
63 ALLOC(atm, atm_t, 1);
64 ALLOC(atm_filt, atm_t, 1);
65 ALLOC(work, double,
66 NP);
67 ALLOC(zs, double,
68 NP);
69
70 /* Read control parameters... */
71 mptrac_read_ctl(argv[1], argc, argv, &ctl);
72 const int ens =
73 (int) scan_ctl(argv[1], argc, argv, "STAT_ENS", -1, "-999", NULL);
74 const double p0 =
75 P(scan_ctl(argv[1], argc, argv, "STAT_Z0", -1, "-1000", NULL));
76 const double p1 =
77 P(scan_ctl(argv[1], argc, argv, "STAT_Z1", -1, "1000", NULL));
78 const double lat0 =
79 scan_ctl(argv[1], argc, argv, "STAT_LAT0", -1, "-1000", NULL);
80 const double lat1 =
81 scan_ctl(argv[1], argc, argv, "STAT_LAT1", -1, "1000", NULL);
82 const double lon0 =
83 scan_ctl(argv[1], argc, argv, "STAT_LON0", -1, "-1000", NULL);
84 const double lon1 =
85 scan_ctl(argv[1], argc, argv, "STAT_LON1", -1, "1000", NULL);
86
87 /* Write info... */
88 LOG(1, "Write air parcel statistics: %s", argv[2]);
89
90 /* Create output file... */
91 if (!(out = fopen(argv[2], "w")))
92 ERRMSG("Cannot create file!");
93
94 /* Write header... */
95 fprintf(out,
96 "# $1 = time [s]\n"
97 "# $2 = time difference [s]\n"
98 "# $3 = altitude (%s) [km]\n"
99 "# $4 = longitude (%s) [deg]\n"
100 "# $5 = latitude (%s) [deg]\n", argv[3], argv[3], argv[3]);
101 for (int iq = 0; iq < ctl.nq; iq++)
102 fprintf(out, "# $%d = %s (%s) [%s]\n", iq + 6,
103 ctl.qnt_name[iq], argv[3], ctl.qnt_unit[iq]);
104 fprintf(out, "# $%d = number of particles\n\n", ctl.nq + 6);
105
106 /* Loop over files... */
107 for (int f = 4; f < argc; f++) {
108
109 /* Read atmopheric data... */
110 if (!mptrac_read_atm(argv[f], &ctl, atm))
111 continue;
112
113 /* Get time from filename... */
114 int time_offset = ctl.atm_type < 2 ? 23 : 22;
115 t = time_from_filename(argv[f], time_offset, 1);
116
117 /* Save initial time... */
118 if (!init) {
119 init = 1;
120 t0 = t;
121 }
122
123 /* Filter data... */
124 atm_filt->np = 0;
125 for (int ip = 0; ip < atm->np; ip++) {
126
127 /* Check time... */
128 if (!isfinite(atm->time[ip]))
129 continue;
130
131 /* Check ensemble index... */
132 if (ctl.qnt_ens > 0 && atm->q[ctl.qnt_ens][ip] != ens)
133 continue;
134
135 /* Check spatial range... */
136 if (atm->p[ip] > p0 || atm->p[ip] < p1
137 || atm->lon[ip] < lon0 || atm->lon[ip] > lon1
138 || atm->lat[ip] < lat0 || atm->lat[ip] > lat1)
139 continue;
140
141 /* Save data... */
142 atm_filt->time[atm_filt->np] = atm->time[ip];
143 atm_filt->p[atm_filt->np] = atm->p[ip];
144 atm_filt->lon[atm_filt->np] = atm->lon[ip];
145 atm_filt->lat[atm_filt->np] = atm->lat[ip];
146 for (int iq = 0; iq < ctl.nq; iq++)
147 atm_filt->q[iq][atm_filt->np] = atm->q[iq][ip];
148 atm_filt->np++;
149 }
150
151 /* Get heights... */
152 for (int ip = 0; ip < atm_filt->np; ip++)
153 zs[ip] = Z(atm_filt->p[ip]);
154
155 /* Get statistics... */
156 if (strcasecmp(argv[3], "mean") == 0) {
157 zm = gsl_stats_mean(zs, 1, (size_t) atm_filt->np);
158 lonm = gsl_stats_mean(atm_filt->lon, 1, (size_t) atm_filt->np);
159 latm = gsl_stats_mean(atm_filt->lat, 1, (size_t) atm_filt->np);
160 for (int iq = 0; iq < ctl.nq; iq++)
161 qm[iq] = gsl_stats_mean(atm_filt->q[iq], 1, (size_t) atm_filt->np);
162 } else if (strcasecmp(argv[3], "stddev") == 0) {
163 zm = gsl_stats_sd(zs, 1, (size_t) atm_filt->np);
164 lonm = gsl_stats_sd(atm_filt->lon, 1, (size_t) atm_filt->np);
165 latm = gsl_stats_sd(atm_filt->lat, 1, (size_t) atm_filt->np);
166 for (int iq = 0; iq < ctl.nq; iq++)
167 qm[iq] = gsl_stats_sd(atm_filt->q[iq], 1, (size_t) atm_filt->np);
168 } else if (strcasecmp(argv[3], "min") == 0) {
169 zm = gsl_stats_min(zs, 1, (size_t) atm_filt->np);
170 lonm = gsl_stats_min(atm_filt->lon, 1, (size_t) atm_filt->np);
171 latm = gsl_stats_min(atm_filt->lat, 1, (size_t) atm_filt->np);
172 for (int iq = 0; iq < ctl.nq; iq++)
173 qm[iq] = gsl_stats_min(atm_filt->q[iq], 1, (size_t) atm_filt->np);
174 } else if (strcasecmp(argv[3], "max") == 0) {
175 zm = gsl_stats_max(zs, 1, (size_t) atm_filt->np);
176 lonm = gsl_stats_max(atm_filt->lon, 1, (size_t) atm_filt->np);
177 latm = gsl_stats_max(atm_filt->lat, 1, (size_t) atm_filt->np);
178 for (int iq = 0; iq < ctl.nq; iq++)
179 qm[iq] = gsl_stats_max(atm_filt->q[iq], 1, (size_t) atm_filt->np);
180 } else if (strcasecmp(argv[3], "skew") == 0) {
181 zm = gsl_stats_skew(zs, 1, (size_t) atm_filt->np);
182 lonm = gsl_stats_skew(atm_filt->lon, 1, (size_t) atm_filt->np);
183 latm = gsl_stats_skew(atm_filt->lat, 1, (size_t) atm_filt->np);
184 for (int iq = 0; iq < ctl.nq; iq++)
185 qm[iq] = gsl_stats_skew(atm_filt->q[iq], 1, (size_t) atm_filt->np);
186 } else if (strcasecmp(argv[3], "kurt") == 0) {
187 zm = gsl_stats_kurtosis(zs, 1, (size_t) atm_filt->np);
188 lonm = gsl_stats_kurtosis(atm_filt->lon, 1, (size_t) atm_filt->np);
189 latm = gsl_stats_kurtosis(atm_filt->lat, 1, (size_t) atm_filt->np);
190 for (int iq = 0; iq < ctl.nq; iq++)
191 qm[iq] =
192 gsl_stats_kurtosis(atm_filt->q[iq], 1, (size_t) atm_filt->np);
193 } else if (strcasecmp(argv[3], "median") == 0) {
194 zm = gsl_stats_median(zs, 1, (size_t) atm_filt->np);
195 lonm = gsl_stats_median(atm_filt->lon, 1, (size_t) atm_filt->np);
196 latm = gsl_stats_median(atm_filt->lat, 1, (size_t) atm_filt->np);
197 for (int iq = 0; iq < ctl.nq; iq++)
198 qm[iq] = gsl_stats_median(atm_filt->q[iq], 1, (size_t) atm_filt->np);
199 } else if (strcasecmp(argv[3], "absdev") == 0) {
200 zm = gsl_stats_absdev(zs, 1, (size_t) atm_filt->np);
201 lonm = gsl_stats_absdev(atm_filt->lon, 1, (size_t) atm_filt->np);
202 latm = gsl_stats_absdev(atm_filt->lat, 1, (size_t) atm_filt->np);
203 for (int iq = 0; iq < ctl.nq; iq++)
204 qm[iq] = gsl_stats_absdev(atm_filt->q[iq], 1, (size_t) atm_filt->np);
205 } else if (strcasecmp(argv[3], "mad") == 0) {
206 zm = gsl_stats_mad0(zs, 1, (size_t) atm_filt->np, work);
207 lonm = gsl_stats_mad0(atm_filt->lon, 1, (size_t) atm_filt->np, work);
208 latm = gsl_stats_mad0(atm_filt->lat, 1, (size_t) atm_filt->np, work);
209 for (int iq = 0; iq < ctl.nq; iq++)
210 qm[iq] =
211 gsl_stats_mad0(atm_filt->q[iq], 1, (size_t) atm_filt->np, work);
212 } else
213 ERRMSG("Unknown parameter!");
214
215 /* Write data... */
216 fprintf(out, "%.2f %.2f %g %g %g", t, t - t0, zm, lonm, latm);
217 for (int iq = 0; iq < ctl.nq; iq++) {
218 fprintf(out, " ");
219 fprintf(out, ctl.qnt_format[iq], qm[iq]);
220 }
221 fprintf(out, " %d\n", atm_filt->np);
222 }
223
224 /* Close file... */
225 fclose(out);
226
227 /* Free... */
228 free(atm);
229 free(atm_filt);
230 free(work);
231 free(zs);
232
233 return EXIT_SUCCESS;
234}
scan_ctl
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.
Definition: mptrac.c:11856
mptrac_read_atm
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.
Definition: mptrac.c:6088
time_from_filename
double time_from_filename(const char *filename, const int offset, const int with_seconds)
Extracts and converts a timestamp from a filename to Julian seconds.
Definition: mptrac.c:12128
mptrac_read_ctl
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.
Definition: mptrac.c:6223
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:2172
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: mptrac.h:1979
Z
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:2009
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1550
NQ
#define NQ
Maximum number of quantities per data point.
Definition: mptrac.h:364
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:457
NP
#define NP
Maximum number of atmospheric data points.
Definition: mptrac.h:359
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:2102
atm_t
Air parcel data.
Definition: mptrac.h:3311
atm_t::time
double time[NP]
Time [s].
Definition: mptrac.h:3317
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3326
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3323
atm_t::np
int np
Number of air parcels.
Definition: mptrac.h:3314
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3329
atm_t::p
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3320
ctl_t
Control parameters.
Definition: mptrac.h:2260
ctl_t::qnt_format
char qnt_format[NQ][LEN]
Quantity output format.
Definition: mptrac.h:2279
ctl_t::atm_type
int atm_type
Type of atmospheric data files (0=ASCII, 1=binary, 2=netCDF, 3=CLaMS_traj, 4=CLaMS_pos).
Definition: mptrac.h:3077
ctl_t::qnt_unit
char qnt_unit[NQ][LEN]
Quantity units.
Definition: mptrac.h:2276
ctl_t::qnt_name
char qnt_name[NQ][LEN]
Quantity names.
Definition: mptrac.h:2270
ctl_t::qnt_ens
int qnt_ens
Quantity array index for ensemble IDs.
Definition: mptrac.h:2285
ctl_t::nq
int nq
Number of quantities.
Definition: mptrac.h:2267
Here is the call graph for this function:
Generated by doxygen 1.9.4