MPTRAC
Loading...
Searching...
No Matches
tropo_sample.c
Go to the documentation of this file.
1/*
2 This file is part of MPTRAC.
3
4 MPTRAC is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 MPTRAC is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
16
17 Copyright (C) 2013-2026 Forschungszentrum Juelich GmbH
18*/
19
25#include "mptrac.h"
26
27/* ------------------------------------------------------------
28 Dimensions...
29 ------------------------------------------------------------ */
30
32#define NT 744
33
34/* ------------------------------------------------------------
35 Functions...
36 ------------------------------------------------------------ */
37
39void usage(
40 void);
41
42/* ------------------------------------------------------------
43 Main...
44 ------------------------------------------------------------ */
45
46int main(
47 int argc,
48 char *argv[]) {
49
50 ctl_t ctl;
51
52 atm_t *atm;
53
54 static FILE *out;
55
56 static char varname[LEN];
57
58 static double times[NT], lons[EX], lats[EY], time0, time1, z0, z0sig, p0,
59 p0sig, t0, t0sig, q0, q0sig, o30, o30sig;
60
61 static float help[EX * EY], tropo_z0[EX][EY], tropo_z1[EX][EY],
62 tropo_p0[EX][EY], tropo_p1[EX][EY], tropo_t0[EX][EY], tropo_t1[EX][EY],
63 tropo_q0[EX][EY], tropo_q1[EX][EY], tropo_o30[EX][EY], tropo_o31[EX][EY];
64
65 static int it_old =
66 -999, ncid, varid, varid_z, varid_p, varid_t, varid_q, varid_o3, h2o, o3,
67 ntime, nlon, nlat, ilon, ilat;
68
69 static size_t count[10], start[10];
70
71 /* Print usage information... */
72 USAGE;
73
74 /* Check arguments... */
75 if (argc < 5)
76 ERRMSG("Missing or invalid command-line arguments.\n\n"
77 "Usage: tropo_sample <ctl> <sample.tab> <tropo.nc> <var> <atm_in>\n\n"
78 "Use -h for full help.");
79
80 /* Allocate... */
81 ALLOC(atm, atm_t, 1);
82
83 /* Read control parameters... */
84 mptrac_read_ctl(argv[1], argc, argv, &ctl);
85 int method =
86 (int) scan_ctl(argv[1], argc, argv, "TROPO_SAMPLE_METHOD", -1, "1", NULL);
87
88 /* Read atmospheric data... */
89 if (!mptrac_read_atm(argv[5], &ctl, atm))
90 ERRMSG("Cannot open file!");
91
92 /* Open tropopause file... */
93 LOG(1, "Read tropopause data: %s", argv[3]);
94 if (nc_open(argv[3], NC_NOWRITE, &ncid) != NC_NOERR)
95 ERRMSG("Cannot open file!");
96
97 /* Get dimensions... */
98 NC_INQ_DIM("time", &ntime, 1, NT, 1);
99 NC_INQ_DIM("lat", &nlat, 1, EY, 1);
100 NC_INQ_DIM("lon", &nlon, 1, EX, 1);
101
102 /* Read coordinates... */
103 NC_GET_DOUBLE("time", times, 1);
104 NC_GET_DOUBLE("lat", lats, 1);
105 NC_GET_DOUBLE("lon", lons, 1);
106
107 /* Get variable indices... */
108 sprintf(varname, "%s_z", argv[4]);
109 NC(nc_inq_varid(ncid, varname, &varid_z));
110 sprintf(varname, "%s_p", argv[4]);
111 NC(nc_inq_varid(ncid, varname, &varid_p));
112 sprintf(varname, "%s_t", argv[4]);
113 NC(nc_inq_varid(ncid, varname, &varid_t));
114 sprintf(varname, "%s_q", argv[4]);
115 h2o = (nc_inq_varid(ncid, varname, &varid_q) == NC_NOERR);
116 sprintf(varname, "%s_o3", argv[4]);
117 o3 = (nc_inq_varid(ncid, varname, &varid_o3) == NC_NOERR);
118
119 /* Set dimensions... */
120 count[0] = 1;
121 count[1] = (size_t) nlat;
122 count[2] = (size_t) nlon;
123
124 /* Create file... */
125 LOG(1, "Write tropopause sample data: %s", argv[2]);
126 if (!(out = fopen(argv[2], "w")))
127 ERRMSG("Cannot create file!");
128
129 /* Write header... */
130 fprintf(out,
131 "# $1 = time [s]\n"
132 "# $2 = altitude [km]\n"
133 "# $3 = longitude [deg]\n" "# $4 = latitude [deg]\n");
134 for (int iq = 0; iq < ctl.nq; iq++)
135 fprintf(out, "# $%i = %s [%s]\n", iq + 5, ctl.qnt_name[iq],
136 ctl.qnt_unit[iq]);
137 fprintf(out, "# $%d = tropopause height (mean) [km]\n", 5 + ctl.nq);
138 fprintf(out, "# $%d = tropopause pressure (mean) [hPa]\n", 6 + ctl.nq);
139 fprintf(out, "# $%d = tropopause temperature (mean) [K]\n", 7 + ctl.nq);
140 fprintf(out, "# $%d = tropopause water vapor (mean) [ppv]\n", 8 + ctl.nq);
141 fprintf(out, "# $%d = tropopause ozone (mean) [ppv]\n", 9 + ctl.nq);
142 fprintf(out, "# $%d = tropopause height (sigma) [km]\n", 10 + ctl.nq);
143 fprintf(out, "# $%d = tropopause pressure (sigma) [hPa]\n", 11 + ctl.nq);
144 fprintf(out, "# $%d = tropopause temperature (sigma) [K]\n", 12 + ctl.nq);
145 fprintf(out, "# $%d = tropopause water vapor (sigma) [ppv]\n", 13 + ctl.nq);
146 fprintf(out, "# $%d = tropopause ozone (sigma) [ppv]\n\n", 14 + ctl.nq);
147
148 /* Loop over particles... */
149 for (int ip = 0; ip < atm->np; ip++) {
150
151 /* Check temporal ordering... */
152 if (ip > 0 && atm->time[ip] < atm->time[ip - 1])
153 ERRMSG("Time must be ascending!");
154
155 /* Check range... */
156 if (atm->time[ip] < times[0] || atm->time[ip] > times[ntime - 1])
157 continue;
158
159 /* Read data... */
160 int it = locate_irr(times, (int) ntime, atm->time[ip]);
161 if (it != it_old) {
162
163 time0 = times[it];
164 start[0] = (size_t) it;
165 NC(nc_get_vara_float(ncid, varid_z, start, count, help));
166 for (ilon = 0; ilon < nlon; ilon++)
167 for (ilat = 0; ilat < nlat; ilat++)
168 tropo_z0[ilon][ilat] = help[ilat * nlon + ilon];
169 NC(nc_get_vara_float(ncid, varid_p, start, count, help));
170 for (ilon = 0; ilon < nlon; ilon++)
171 for (ilat = 0; ilat < nlat; ilat++)
172 tropo_p0[ilon][ilat] = help[ilat * nlon + ilon];
173 NC(nc_get_vara_float(ncid, varid_t, start, count, help));
174 for (ilon = 0; ilon < nlon; ilon++)
175 for (ilat = 0; ilat < nlat; ilat++)
176 tropo_t0[ilon][ilat] = help[ilat * nlon + ilon];
177 if (h2o) {
178 NC(nc_get_vara_float(ncid, varid_q, start, count, help));
179 for (ilon = 0; ilon < nlon; ilon++)
180 for (ilat = 0; ilat < nlat; ilat++)
181 tropo_q0[ilon][ilat] = help[ilat * nlon + ilon];
182 } else
183 for (ilon = 0; ilon < nlon; ilon++)
184 for (ilat = 0; ilat < nlat; ilat++)
185 tropo_q0[ilon][ilat] = NAN;
186 if (o3) {
187 NC(nc_get_vara_float(ncid, varid_o3, start, count, help));
188 for (ilon = 0; ilon < nlon; ilon++)
189 for (ilat = 0; ilat < nlat; ilat++)
190 tropo_o30[ilon][ilat] = help[ilat * nlon + ilon];
191 } else
192 for (ilon = 0; ilon < nlon; ilon++)
193 for (ilat = 0; ilat < nlat; ilat++)
194 tropo_o30[ilon][ilat] = NAN;
195
196 time1 = times[it + 1];
197 start[0] = (size_t) it + 1;
198 NC(nc_get_vara_float(ncid, varid_z, start, count, help));
199 for (ilon = 0; ilon < nlon; ilon++)
200 for (ilat = 0; ilat < nlat; ilat++)
201 tropo_z1[ilon][ilat] = help[ilat * nlon + ilon];
202 NC(nc_get_vara_float(ncid, varid_p, start, count, help));
203 for (ilon = 0; ilon < nlon; ilon++)
204 for (ilat = 0; ilat < nlat; ilat++)
205 tropo_p1[ilon][ilat] = help[ilat * nlon + ilon];
206 NC(nc_get_vara_float(ncid, varid_t, start, count, help));
207 for (ilon = 0; ilon < nlon; ilon++)
208 for (ilat = 0; ilat < nlat; ilat++)
209 tropo_t1[ilon][ilat] = help[ilat * nlon + ilon];
210 if (h2o) {
211 NC(nc_get_vara_float(ncid, varid_q, start, count, help));
212 for (ilon = 0; ilon < nlon; ilon++)
213 for (ilat = 0; ilat < nlat; ilat++)
214 tropo_q1[ilon][ilat] = help[ilat * nlon + ilon];
215 } else
216 for (ilon = 0; ilon < nlon; ilon++)
217 for (ilat = 0; ilat < nlat; ilat++)
218 tropo_q1[ilon][ilat] = NAN;
219 if (o3) {
220 NC(nc_get_vara_float(ncid, varid_o3, start, count, help));
221 for (ilon = 0; ilon < nlon; ilon++)
222 for (ilat = 0; ilat < nlat; ilat++)
223 tropo_o31[ilon][ilat] = help[ilat * nlon + ilon];
224 } else
225 for (ilon = 0; ilon < nlon; ilon++)
226 for (ilat = 0; ilat < nlat; ilat++)
227 tropo_o31[ilon][ilat] = NAN;
228 }
229 it_old = it;
230
231 /* Interpolate... */
232 intpol_tropo_3d(time0, tropo_z0, time1, tropo_z1,
233 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
234 atm->lat[ip], method, &z0, &z0sig);
235 intpol_tropo_3d(time0, tropo_p0, time1, tropo_p1,
236 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
237 atm->lat[ip], method, &p0, &p0sig);
238 intpol_tropo_3d(time0, tropo_t0, time1, tropo_t1,
239 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
240 atm->lat[ip], method, &t0, &t0sig);
241 intpol_tropo_3d(time0, tropo_q0, time1, tropo_q1,
242 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
243 atm->lat[ip], method, &q0, &q0sig);
244 intpol_tropo_3d(time0, tropo_o30, time1, tropo_o31,
245 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
246 atm->lat[ip], method, &o30, &o30sig);
247
248 /* Write output... */
249 fprintf(out, "%.2f %g %g %g", atm->time[ip], Z(atm->p[ip]),
250 atm->lon[ip], atm->lat[ip]);
251 for (int iq = 0; iq < ctl.nq; iq++) {
252 fprintf(out, " ");
253 fprintf(out, ctl.qnt_format[iq], atm->q[iq][ip]);
254 }
255 fprintf(out, " %g %g %g %g %g %g %g %g %g %g\n",
256 z0, p0, t0, q0, o30, z0sig, p0sig, t0sig, q0sig, o30sig);
257 }
258
259 /* Close files... */
260 fclose(out);
261 NC(nc_close(ncid));
262
263 /* Free... */
264 free(atm);
265
266 return EXIT_SUCCESS;
267}
268
269/*****************************************************************************/
270
272void usage(
273 void) {
274
275 printf("\nMPTRAC tropo_sample tool.\n\n");
276 printf("Sample tropopause data at atmospheric particle locations.\n");
277 printf("\n");
278 printf("Usage:\n");
279 printf(" tropo_sample <ctl> <sample.tab> <tropo.nc> <var> <atm_in>\n");
280 printf("\n");
281 printf("Arguments:\n");
282 printf(" <ctl> Control file.\n");
283 printf(" <sample.tab> Output table.\n");
284 printf(" <tropo.nc> Tropopause netCDF input file.\n");
285 printf
286 (" <var> Tropopause variable prefix, such as clp, dyn, wmo_1st,\n");
287 printf(" or wmo_2nd.\n");
288 printf(" <atm_in> Atmospheric input file with sample locations.\n");
289 printf("\nFurther information:\n");
290 printf(" Manual: https://slcs-jsc.github.io/mptrac/\n");
291}
locate_irr
int locate_irr(const double *xx, const int n, const double x)
Locate the index of the interval containing a given value in a sorted array.
Definition mptrac.c:3801
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:13329
intpol_tropo_3d
void intpol_tropo_3d(const double time0, float array0[EX][EY], const double time1, float array1[EX][EY], const double lons[EX], const double lats[EY], const int nlon, const int nlat, const double time, const double lon, const double lat, const int method, double *var, double *sigma)
Interpolates tropopause data in 3D (latitude, longitude, and time).
Definition mptrac.c:3480
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:6954
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:7089
mptrac.h
MPTRAC library declarations.
LEN
#define LEN
Maximum length of ASCII data lines.
Definition mptrac.h:559
NC
#define NC(cmd)
Execute a NetCDF command and check for errors.
Definition mptrac.h:1517
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition mptrac.h:2426
EY
#define EY
Maximum number of latitudes for meteo data.
Definition mptrac.h:554
USAGE
#define USAGE
Print usage information on -h or --help.
Definition mptrac.h:2233
Z
#define Z(p)
Convert pressure to altitude.
Definition mptrac.h:2263
EX
#define EX
Maximum number of longitudes for meteo data.
Definition mptrac.h:549
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition mptrac.h:668
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition mptrac.h:2356
NC_GET_DOUBLE
#define NC_GET_DOUBLE(varname, ptr, force)
Retrieve a double-precision variable from a NetCDF file.
Definition mptrac.h:1587
NC_INQ_DIM
#define NC_INQ_DIM(dimname, ptr, min, max, check)
Inquire the length of a dimension in a NetCDF file.
Definition mptrac.h:1617
atm_t
Air parcel data.
Definition mptrac.h:3612
atm_t::time
double time[NP]
Time [s].
Definition mptrac.h:3618
atm_t::lat
double lat[NP]
Latitude [deg].
Definition mptrac.h:3627
atm_t::lon
double lon[NP]
Longitude [deg].
Definition mptrac.h:3624
atm_t::np
int np
Number of air parcels.
Definition mptrac.h:3615
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition mptrac.h:3630
atm_t::p
double p[NP]
Pressure [hPa].
Definition mptrac.h:3621
ctl_t
Control parameters.
Definition mptrac.h:2514
ctl_t::qnt_format
char qnt_format[NQ][LEN]
Quantity output format.
Definition mptrac.h:2533
ctl_t::qnt_unit
char qnt_unit[NQ][LEN]
Quantity units.
Definition mptrac.h:2530
ctl_t::qnt_name
char qnt_name[NQ][LEN]
Quantity names.
Definition mptrac.h:2524
ctl_t::nq
int nq
Number of quantities.
Definition mptrac.h:2521
main
int main(int argc, char *argv[])
Definition tropo_sample.c:46
NT
#define NT
Maximum number of time steps.
Definition tropo_sample.c:32
usage
void usage(void)
Print command-line help.
Definition tropo_sample.c:272
Generated on Tue Nov 10 2026 15:48:54 for MPTRAC by doxygen 1.12.0