MPTRAC
Functions
atm_select.c File Reference

Extract subsets of air parcels from atmospheric data files. More...

#include "mptrac.h"

Go to the source code of this file.

Functions

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

Detailed Description

Extract subsets of air parcels from atmospheric data files.

Definition in file atm_select.c.

Function Documentation

◆ main()

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

Definition at line 27 of file atm_select.c.

29 {
30
31 ctl_t ctl;
32
33 atm_t *atm, *atm2;
34
35 /* Allocate... */
36 ALLOC(atm, atm_t, 1);
37 ALLOC(atm2, atm_t, 1);
38
39 /* Check arguments... */
40 if (argc < 4)
41 ERRMSG("Give parameters: <ctl> <atm_select> <atm1> [<atm2> ...]");
42
43 /* Read control parameters... */
44 read_ctl(argv[1], argc, argv, &ctl);
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 /* Get Cartesian coordinates... */
73 double x0[3], x1[3];
74 geo2cart(0, rlon, rlat, x0);
75
76 /* Loop over files... */
77 for (int f = 3; f < argc; f++) {
78
79 /* Read atmopheric data... */
80 if (!read_atm(argv[f], &ctl, atm))
81 continue;
82
83 /* Adjust range of air parcels... */
84 if (ip0 < 0)
85 ip0 = 0;
86 ip0 = MIN(ip0, atm->np - 1);
87 if (ip1 < 0)
88 ip1 = atm->np - 1;
89 ip1 = MIN(ip1, atm->np - 1);
90 if (ip1 < ip0)
91 ip1 = ip0;
92
93 /* Loop over air parcels... */
94 for (int ip = ip0; ip <= ip1; ip += stride) {
95
96 /* Check air parcel index... */
97 if (ctl.qnt_idx >= 0 && idx0 >= 0 && idx1 >= 0)
98 if (atm->q[ctl.qnt_idx][ip] < idx0 || atm->q[ctl.qnt_idx][ip] > idx1)
99 continue;
100
101 /* Check time... */
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 /* Check vertical distance... */
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 /* Check potential temperature... */
114 if (theta0 != theta1) {
115 double theta;
116 if (ctl.qnt_theta >= 0)
117 theta = atm->q[ctl.qnt_theta][ip];
118 else if (ctl.qnt_t >= 0)
119 theta = THETA(atm->p[ip], atm->q[ctl.qnt_t][ip]);
120 else
121 ERRMSG
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 /* Check longitude... */
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 /* Check latitude... */
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 /* Check horizontal distace... */
141 if (r0 != r1) {
142 geo2cart(0, atm->lon[ip], atm->lat[ip], x1);
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 /* Copy data... */
150 atm2->time[atm2->np] = atm->time[ip];
151 atm2->p[atm2->np] = atm->p[ip];
152 atm2->lon[atm2->np] = atm->lon[ip];
153 atm2->lat[atm2->np] = atm->lat[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 /* Close file... */
162 write_atm(argv[2], &ctl, atm2, 0);
163
164 /* Free... */
165 free(atm);
166 free(atm2);
167
168 return EXIT_SUCCESS;
169}
int read_atm(const char *filename, ctl_t *ctl, atm_t *atm)
Reads air parcel data from a specified file into the given atmospheric structure.
Definition: mptrac.c:4183
void write_atm(const char *filename, ctl_t *ctl, atm_t *atm, double t)
Writes air parcel data to a file in various formats.
Definition: mptrac.c:8317
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.
Definition: mptrac.c:4789
double scan_ctl(const char *filename, int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Scans a control file or command-line arguments for a specified variable.
Definition: mptrac.c:7943
void geo2cart(const double z, const double lon, const double lat, double *x)
Converts geographic coordinates (longitude, latitude, altitude) to Cartesian coordinates.
Definition: mptrac.c:941
#define MIN(a, b)
Macro to determine the minimum of two values.
Definition: mptrac.h:952
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:1881
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1255
#define THETA(p, t)
Compute potential temperature.
Definition: mptrac.h:1609
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:344
#define NP
Maximum number of atmospheric data points.
Definition: mptrac.h:241
#define DIST(a, b)
Calculate the distance between two points in Cartesian coordinates.
Definition: mptrac.h:552
Air parcel data.
Definition: mptrac.h:3068
double time[NP]
Time [s].
Definition: mptrac.h:3074
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3083
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3080
int np
Number of air parcels.
Definition: mptrac.h:3071
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3086
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3077
Control parameters.
Definition: mptrac.h:2135
int qnt_theta
Quantity array index for potential temperature.
Definition: mptrac.h:2355
int qnt_t
Quantity array index for temperature.
Definition: mptrac.h:2238
int qnt_idx
Quantity array index for air parcel IDs.
Definition: mptrac.h:2175
int nq
Number of quantities.
Definition: mptrac.h:2160
Here is the call graph for this function: