29 {
30
32
34
35
38
39
40 if (argc < 4)
41 ERRMSG(
"Give parameters: <ctl> <atm_select> <atm1> [<atm2> ...]");
42
43
45 int stride =
46 (int)
scan_ctl(argv[1], argc, argv,
"SELECT_STRIDE", -1,
"1", NULL);
47 int idx0 =
48 (int)
scan_ctl(argv[1], argc, argv,
"SELECT_IDX0", -1,
"-999", NULL);
49 int idx1 =
50 (int)
scan_ctl(argv[1], argc, argv,
"SELECT_IDX1", -1,
"-999", NULL);
51 int ip0 =
52 (int)
scan_ctl(argv[1], argc, argv,
"SELECT_IP0", -1,
"-999", NULL);
53 int ip1 =
54 (int)
scan_ctl(argv[1], argc, argv,
"SELECT_IP1", -1,
"-999", NULL);
55 double t0 =
scan_ctl(argv[1], argc, argv,
"SELECT_T0", -1,
"0", NULL);
56 double t1 =
scan_ctl(argv[1], argc, argv,
"SELECT_T1", -1,
"0", NULL);
57 double p0 =
P(
scan_ctl(argv[1], argc, argv,
"SELECT_Z0", -1,
"0", NULL));
58 double p1 =
P(
scan_ctl(argv[1], argc, argv,
"SELECT_Z1", -1,
"0", NULL));
59 double theta0 =
60 scan_ctl(argv[1], argc, argv,
"SELECT_THETA0", -1,
"0", NULL);
61 double theta1 =
62 scan_ctl(argv[1], argc, argv,
"SELECT_THETA1", -1,
"0", NULL);
63 double lon0 =
scan_ctl(argv[1], argc, argv,
"SELECT_LON0", -1,
"0", NULL);
64 double lon1 =
scan_ctl(argv[1], argc, argv,
"SELECT_LON1", -1,
"0", NULL);
65 double lat0 =
scan_ctl(argv[1], argc, argv,
"SELECT_LAT0", -1,
"0", NULL);
66 double lat1 =
scan_ctl(argv[1], argc, argv,
"SELECT_LAT1", -1,
"0", NULL);
67 double r0 =
scan_ctl(argv[1], argc, argv,
"SELECT_R0", -1,
"0", NULL);
68 double r1 =
scan_ctl(argv[1], argc, argv,
"SELECT_R1", -1,
"0", NULL);
69 double rlon =
scan_ctl(argv[1], argc, argv,
"SELECT_RLON", -1,
"0", NULL);
70 double rlat =
scan_ctl(argv[1], argc, argv,
"SELECT_RLAT", -1,
"0", NULL);
71
72
73 double x0[3], x1[3];
75
76
77 for (int f = 3; f < argc; f++) {
78
79
81 continue;
82
83
84 if (ip0 < 0)
85 ip0 = 0;
86 ip0 =
MIN(ip0, atm->
np - 1);
87 if (ip1 < 0)
89 ip1 =
MIN(ip1, atm->
np - 1);
90 if (ip1 < ip0)
91 ip1 = ip0;
92
93
94 for (int ip = ip0; ip <= ip1; ip += stride) {
95
96
97 if (ctl.
qnt_idx >= 0 && idx0 >= 0 && idx1 >= 0)
99 continue;
100
101
102 if (t0 != t1)
103 if ((t1 > t0 && (atm->
time[ip] < t0 || atm->
time[ip] > t1))
104 || (t1 < t0 && (atm->
time[ip] < t0 && atm->
time[ip] > t1)))
105 continue;
106
107
108 if (p0 != p1)
109 if ((p0 > p1 && (atm->
p[ip] > p0 || atm->
p[ip] < p1))
110 || (p0 < p1 && (atm->
p[ip] > p0 && atm->
p[ip] < p1)))
111 continue;
112
113
114 if (theta0 != theta1) {
115 double theta;
118 else if (ctl.
qnt_t >= 0)
120 else
122 ("Filtering requires temperature or potential temperature data!");
123 if ((theta1 > theta0 && (theta < theta0 || theta > theta1))
124 || (theta1 < theta0 && (theta < theta0 && theta > theta1)))
125 continue;
126 }
127
128
129 if (lon0 != lon1)
130 if ((lon1 > lon0 && (atm->
lon[ip] < lon0 || atm->
lon[ip] > lon1))
131 || (lon1 < lon0 && (atm->
lon[ip] < lon0 && atm->
lon[ip] > lon1)))
132 continue;
133
134
135 if (lat0 != lat1)
136 if ((lat1 > lat0 && (atm->
lat[ip] < lat0 || atm->
lat[ip] > lat1))
137 || (lat1 < lat0 && (atm->
lat[ip] < lat0 && atm->
lat[ip] > lat1)))
138 continue;
139
140
141 if (r0 != r1) {
143 double r =
DIST(x0, x1);
144 if ((r1 > r0 && (r < r0 || r > r1))
145 || (r1 < r0 && (r < r0 && r > r1)))
146 continue;
147 }
148
149
151 atm2->
p[atm2->
np] = atm->
p[ip];
154 for (
int iq = 0; iq < ctl.
nq; iq++)
155 atm2->
q[iq][atm2->
np] = atm->
q[iq][ip];
156 if ((++atm2->
np) >
NP)
157 ERRMSG(
"Too many air parcels!");
158 }
159 }
160
161
163
164
165 free(atm);
166 free(atm2);
167
168 return EXIT_SUCCESS;
169}
void write_atm(const char *filename, const ctl_t *ctl, const atm_t *atm, const double t)
Writes air parcel data to a file in various formats.
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.
void read_ctl(const char *filename, int argc, char *argv[], ctl_t *ctl)
Reads control parameters from a configuration file and populates the given structure.
int 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 geo2cart(const double z, const double lon, const double lat, double *x)
Converts geographic coordinates (longitude, latitude, altitude) to Cartesian coordinates.
#define MIN(a, b)
Macro to determine the minimum of two values.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
#define P(z)
Compute pressure at given altitude.
#define THETA(p, t)
Compute potential temperature.
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
#define NP
Maximum number of atmospheric data points.
#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].
int qnt_theta
Quantity array index for potential temperature.
int qnt_t
Quantity array index for temperature.
int qnt_idx
Quantity array index for air parcel IDs.
int nq
Number of quantities.