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 238 of file atm_stat.c.

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

◆ 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 /* Allocate... */
54 ALLOC(atm, atm_t, 1);
55 ALLOC(atm_filt, atm_t, 1);
56 ALLOC(work, double,
57 NP);
58 ALLOC(zs, double,
59 NP);
60
61 /* Print usage information... */
62 USAGE;
63
64 /* Check arguments... */
65 if (argc < 4)
66 ERRMSG("Missing or invalid command-line arguments.\n\n"
67 "Usage: atm_stat <ctl> <stat.tab> <param> <atm1> [<atm2> ...]\n\n"
68 "Use -h for full help.");
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 t = time_from_filename(argv[f], ctl.atm_type < 2 ? 20 : 19);
115
116 /* Save initial time... */
117 if (!init) {
118 init = 1;
119 t0 = t;
120 }
121
122 /* Filter data... */
123 atm_filt->np = 0;
124 for (int ip = 0; ip < atm->np; ip++) {
125
126 /* Check time... */
127 if (!isfinite(atm->time[ip]))
128 continue;
129
130 /* Check ensemble index... */
131 if (ctl.qnt_ens > 0 && atm->q[ctl.qnt_ens][ip] != ens)
132 continue;
133
134 /* Check spatial range... */
135 if (atm->p[ip] > p0 || atm->p[ip] < p1
136 || atm->lon[ip] < lon0 || atm->lon[ip] > lon1
137 || atm->lat[ip] < lat0 || atm->lat[ip] > lat1)
138 continue;
139
140 /* Save data... */
141 atm_filt->time[atm_filt->np] = atm->time[ip];
142 atm_filt->p[atm_filt->np] = atm->p[ip];
143 atm_filt->lon[atm_filt->np] = atm->lon[ip];
144 atm_filt->lat[atm_filt->np] = atm->lat[ip];
145 for (int iq = 0; iq < ctl.nq; iq++)
146 atm_filt->q[iq][atm_filt->np] = atm->q[iq][ip];
147 atm_filt->np++;
148 }
149
150 /* Get heights... */
151 for (int ip = 0; ip < atm_filt->np; ip++)
152 zs[ip] = Z(atm_filt->p[ip]);
153
154 /* Get statistics... */
155 if (strcasecmp(argv[3], "mean") == 0) {
156 zm = gsl_stats_mean(zs, 1, (size_t) atm_filt->np);
157 lonm = gsl_stats_mean(atm_filt->lon, 1, (size_t) atm_filt->np);
158 latm = gsl_stats_mean(atm_filt->lat, 1, (size_t) atm_filt->np);
159 for (int iq = 0; iq < ctl.nq; iq++)
160 qm[iq] = gsl_stats_mean(atm_filt->q[iq], 1, (size_t) atm_filt->np);
161 } else if (strcasecmp(argv[3], "stddev") == 0) {
162 zm = gsl_stats_sd(zs, 1, (size_t) atm_filt->np);
163 lonm = gsl_stats_sd(atm_filt->lon, 1, (size_t) atm_filt->np);
164 latm = gsl_stats_sd(atm_filt->lat, 1, (size_t) atm_filt->np);
165 for (int iq = 0; iq < ctl.nq; iq++)
166 qm[iq] = gsl_stats_sd(atm_filt->q[iq], 1, (size_t) atm_filt->np);
167 } else if (strcasecmp(argv[3], "min") == 0) {
168 zm = gsl_stats_min(zs, 1, (size_t) atm_filt->np);
169 lonm = gsl_stats_min(atm_filt->lon, 1, (size_t) atm_filt->np);
170 latm = gsl_stats_min(atm_filt->lat, 1, (size_t) atm_filt->np);
171 for (int iq = 0; iq < ctl.nq; iq++)
172 qm[iq] = gsl_stats_min(atm_filt->q[iq], 1, (size_t) atm_filt->np);
173 } else if (strcasecmp(argv[3], "max") == 0) {
174 zm = gsl_stats_max(zs, 1, (size_t) atm_filt->np);
175 lonm = gsl_stats_max(atm_filt->lon, 1, (size_t) atm_filt->np);
176 latm = gsl_stats_max(atm_filt->lat, 1, (size_t) atm_filt->np);
177 for (int iq = 0; iq < ctl.nq; iq++)
178 qm[iq] = gsl_stats_max(atm_filt->q[iq], 1, (size_t) atm_filt->np);
179 } else if (strcasecmp(argv[3], "skew") == 0) {
180 zm = gsl_stats_skew(zs, 1, (size_t) atm_filt->np);
181 lonm = gsl_stats_skew(atm_filt->lon, 1, (size_t) atm_filt->np);
182 latm = gsl_stats_skew(atm_filt->lat, 1, (size_t) atm_filt->np);
183 for (int iq = 0; iq < ctl.nq; iq++)
184 qm[iq] = gsl_stats_skew(atm_filt->q[iq], 1, (size_t) atm_filt->np);
185 } else if (strcasecmp(argv[3], "kurt") == 0) {
186 zm = gsl_stats_kurtosis(zs, 1, (size_t) atm_filt->np);
187 lonm = gsl_stats_kurtosis(atm_filt->lon, 1, (size_t) atm_filt->np);
188 latm = gsl_stats_kurtosis(atm_filt->lat, 1, (size_t) atm_filt->np);
189 for (int iq = 0; iq < ctl.nq; iq++)
190 qm[iq] =
191 gsl_stats_kurtosis(atm_filt->q[iq], 1, (size_t) atm_filt->np);
192 } else if (strcasecmp(argv[3], "median") == 0) {
193 zm = gsl_stats_median(zs, 1, (size_t) atm_filt->np);
194 lonm = gsl_stats_median(atm_filt->lon, 1, (size_t) atm_filt->np);
195 latm = gsl_stats_median(atm_filt->lat, 1, (size_t) atm_filt->np);
196 for (int iq = 0; iq < ctl.nq; iq++)
197 qm[iq] = gsl_stats_median(atm_filt->q[iq], 1, (size_t) atm_filt->np);
198 } else if (strcasecmp(argv[3], "absdev") == 0) {
199 zm = gsl_stats_absdev(zs, 1, (size_t) atm_filt->np);
200 lonm = gsl_stats_absdev(atm_filt->lon, 1, (size_t) atm_filt->np);
201 latm = gsl_stats_absdev(atm_filt->lat, 1, (size_t) atm_filt->np);
202 for (int iq = 0; iq < ctl.nq; iq++)
203 qm[iq] = gsl_stats_absdev(atm_filt->q[iq], 1, (size_t) atm_filt->np);
204 } else if (strcasecmp(argv[3], "mad") == 0) {
205 zm = gsl_stats_mad0(zs, 1, (size_t) atm_filt->np, work);
206 lonm = gsl_stats_mad0(atm_filt->lon, 1, (size_t) atm_filt->np, work);
207 latm = gsl_stats_mad0(atm_filt->lat, 1, (size_t) atm_filt->np, work);
208 for (int iq = 0; iq < ctl.nq; iq++)
209 qm[iq] =
210 gsl_stats_mad0(atm_filt->q[iq], 1, (size_t) atm_filt->np, work);
211 } else
212 ERRMSG("Unknown parameter!");
213
214 /* Write data... */
215 fprintf(out, "%.2f %.2f %g %g %g", t, t - t0, zm, lonm, latm);
216 for (int iq = 0; iq < ctl.nq; iq++) {
217 fprintf(out, " ");
218 fprintf(out, ctl.qnt_format[iq], qm[iq]);
219 }
220 fprintf(out, " %d\n", atm_filt->np);
221 }
222
223 /* Close file... */
224 fclose(out);
225
226 /* Free... */
227 free(atm);
228 free(atm_filt);
229 free(work);
230 free(zs);
231
232 return EXIT_SUCCESS;
233}
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:10745
time_from_filename
double time_from_filename(const char *filename, const int offset)
Extracts and converts a timestamp from a filename to Julian seconds.
Definition: mptrac.c:11017
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:5117
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:5248
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:2102
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: mptrac.h:1909
Z
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:1939
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1480
NQ
#define NQ
Maximum number of quantities per data point.
Definition: mptrac.h:360
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:453
NP
#define NP
Maximum number of atmospheric data points.
Definition: mptrac.h:355
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:2032
atm_t
Air parcel data.
Definition: mptrac.h:3241
atm_t::time
double time[NP]
Time [s].
Definition: mptrac.h:3247
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3256
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3253
atm_t::np
int np
Number of air parcels.
Definition: mptrac.h:3244
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3259
atm_t::p
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3250
ctl_t
Control parameters.
Definition: mptrac.h:2190
ctl_t::qnt_format
char qnt_format[NQ][LEN]
Quantity output format.
Definition: mptrac.h:2209
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:3007
ctl_t::qnt_unit
char qnt_unit[NQ][LEN]
Quantity units.
Definition: mptrac.h:2206
ctl_t::qnt_name
char qnt_name[NQ][LEN]
Quantity names.
Definition: mptrac.h:2200
ctl_t::qnt_ens
int qnt_ens
Quantity array index for ensemble IDs.
Definition: mptrac.h:2215
ctl_t::nq
int nq
Number of quantities.
Definition: mptrac.h:2197
Here is the call graph for this function:
Generated by doxygen 1.9.4