MPTRAC
tropo_zm.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 static FILE *out;
53
54 static char varname[LEN];
55
56 static double time0, lons[EX], lats[EY], zm[EY], zs[EY], pm[EY],
57 ps[EY], tm[EY], ts[EY], qm[EY], qs[EY], o3m[EY], o3s[EY];
58
59 static float help[EX * EY], tropo_z0[EX][EY], tropo_p0[EX][EY],
60 tropo_t0[EX][EY], tropo_q0[EX][EY], tropo_o30[EX][EY];
61
62 static int ncid, varid, varid_z, varid_p, varid_t, varid_q, varid_o3, h2o,
63 o3, n[EY], nt[EY], init, ntime, nlon, nlat, ilon, ilat;
64
65 static size_t count[10], start[10];
66
67 /* Print usage information... */
68 USAGE;
69
70 /* Check arguments... */
71 if (argc < 5)
72 ERRMSG("Missing or invalid command-line arguments.\n\n"
73 "Usage: tropo_zm <ctl> <zm.tab> <var> <tropo.nc>\n\n"
74 "Use -h for full help.");
75
76 /* Read control parameters... */
77 mptrac_read_ctl(argv[1], argc, argv, &ctl);
78
79 /* Loop over tropopause files... */
80 for (int iarg = 4; iarg < argc; iarg++) {
81
82 /* Open tropopause file... */
83 LOG(1, "Read tropopause data: %s", argv[iarg]);
84 if (nc_open(argv[iarg], NC_NOWRITE, &ncid) != NC_NOERR)
85 ERRMSG("Cannot open file!");
86
87 /* Get dimensions... */
88 NC_INQ_DIM("time", &ntime, 1, NT, 1);
89 NC_INQ_DIM("lat", &nlat, 1, EY, 1);
90 NC_INQ_DIM("lon", &nlon, 1, EX, 1);
91
92 /* Read coordinates... */
93 NC_GET_DOUBLE("lat", lats, 1);
94 NC_GET_DOUBLE("lon", lons, 1);
95
96 /* Get variable indices... */
97 sprintf(varname, "%s_z", argv[3]);
98 NC(nc_inq_varid(ncid, varname, &varid_z));
99 sprintf(varname, "%s_p", argv[3]);
100 NC(nc_inq_varid(ncid, varname, &varid_p));
101 sprintf(varname, "%s_t", argv[3]);
102 NC(nc_inq_varid(ncid, varname, &varid_t));
103 sprintf(varname, "%s_q", argv[3]);
104 h2o = (nc_inq_varid(ncid, varname, &varid_q) == NC_NOERR);
105 sprintf(varname, "%s_o3", argv[3]);
106 o3 = (nc_inq_varid(ncid, varname, &varid_o3) == NC_NOERR);
107
108 /* Set dimensions... */
109 count[0] = 1;
110 count[1] = (size_t) nlat;
111 count[2] = (size_t) nlon;
112
113 /* Loop over time steps... */
114 for (int it = 0; it < ntime; it++) {
115
116 /* Get time from filename... */
117 if (!init) {
118 init = 1;
119 time0 = time_from_filename(argv[iarg], 13);
120 }
121
122 /* Read data... */
123 start[0] = (size_t) it;
124 NC(nc_get_vara_float(ncid, varid_z, start, count, help));
125 for (ilon = 0; ilon < nlon; ilon++)
126 for (ilat = 0; ilat < nlat; ilat++)
127 tropo_z0[ilon][ilat] = help[ilat * nlon + ilon];
128 NC(nc_get_vara_float(ncid, varid_p, start, count, help));
129 for (ilon = 0; ilon < nlon; ilon++)
130 for (ilat = 0; ilat < nlat; ilat++)
131 tropo_p0[ilon][ilat] = help[ilat * nlon + ilon];
132 NC(nc_get_vara_float(ncid, varid_t, start, count, help));
133 for (ilon = 0; ilon < nlon; ilon++)
134 for (ilat = 0; ilat < nlat; ilat++)
135 tropo_t0[ilon][ilat] = help[ilat * nlon + ilon];
136 if (h2o) {
137 NC(nc_get_vara_float(ncid, varid_q, start, count, help));
138 for (ilon = 0; ilon < nlon; ilon++)
139 for (ilat = 0; ilat < nlat; ilat++)
140 tropo_q0[ilon][ilat] = help[ilat * nlon + ilon];
141 } else
142 for (ilon = 0; ilon < nlon; ilon++)
143 for (ilat = 0; ilat < nlat; ilat++)
144 tropo_q0[ilon][ilat] = NAN;
145 if (o3) {
146 NC(nc_get_vara_float(ncid, varid_o3, start, count, help));
147 for (ilon = 0; ilon < nlon; ilon++)
148 for (ilat = 0; ilat < nlat; ilat++)
149 tropo_o30[ilon][ilat] = help[ilat * nlon + ilon];
150 } else
151 for (ilon = 0; ilon < nlon; ilon++)
152 for (ilat = 0; ilat < nlat; ilat++)
153 tropo_o30[ilon][ilat] = NAN;
154
155 /* Averaging... */
156 for (ilat = 0; ilat < nlat; ilat++)
157 for (ilon = 0; ilon < nlon; ilon++) {
158 nt[ilat]++;
159 if (isfinite(tropo_z0[ilon][ilat])
160 && isfinite(tropo_p0[ilon][ilat])
161 && isfinite(tropo_t0[ilon][ilat])
162 && (!h2o || isfinite(tropo_q0[ilon][ilat]))
163 && (!o3 || isfinite(tropo_o30[ilon][ilat]))) {
164 zm[ilat] += tropo_z0[ilon][ilat];
165 zs[ilat] += SQR(tropo_z0[ilon][ilat]);
166 pm[ilat] += tropo_p0[ilon][ilat];
167 ps[ilat] += SQR(tropo_p0[ilon][ilat]);
168 tm[ilat] += tropo_t0[ilon][ilat];
169 ts[ilat] += SQR(tropo_t0[ilon][ilat]);
170 qm[ilat] += tropo_q0[ilon][ilat];
171 qs[ilat] += SQR(tropo_q0[ilon][ilat]);
172 o3m[ilat] += tropo_o30[ilon][ilat];
173 o3s[ilat] += SQR(tropo_o30[ilon][ilat]);
174 n[ilat]++;
175 }
176 }
177 }
178
179 /* Close files... */
180 NC(nc_close(ncid));
181 }
182
183 /* Normalize... */
184 for (ilat = 0; ilat < nlat; ilat++)
185 if (n[ilat] > 0) {
186 zm[ilat] /= n[ilat];
187 pm[ilat] /= n[ilat];
188 tm[ilat] /= n[ilat];
189 qm[ilat] /= n[ilat];
190 o3m[ilat] /= n[ilat];
191 double aux = zs[ilat] / n[ilat] - SQR(zm[ilat]);
192 zs[ilat] = aux > 0 ? sqrt(aux) : 0.0;
193 aux = ps[ilat] / n[ilat] - SQR(pm[ilat]);
194 ps[ilat] = aux > 0 ? sqrt(aux) : 0.0;
195 aux = ts[ilat] / n[ilat] - SQR(tm[ilat]);
196 ts[ilat] = aux > 0 ? sqrt(aux) : 0.0;
197 aux = qs[ilat] / n[ilat] - SQR(qm[ilat]);
198 qs[ilat] = aux > 0 ? sqrt(aux) : 0.0;
199 aux = o3s[ilat] / n[ilat] - SQR(o3m[ilat]);
200 o3s[ilat] = aux > 0 ? sqrt(aux) : 0.0;
201 }
202
203 /* Create file... */
204 LOG(1, "Write tropopause zonal mean data: %s", argv[2]);
205 if (!(out = fopen(argv[2], "w")))
206 ERRMSG("Cannot create file!");
207
208 /* Write header... */
209 fprintf(out,
210 "# $1 = time [s]\n"
211 "# $2 = latitude [deg]\n"
212 "# $3 = tropopause height (mean) [km]\n"
213 "# $4 = tropopause pressure (mean) [hPa]\n"
214 "# $5 = tropopause temperature (mean) [K]\n"
215 "# $6 = tropopause water vapor (mean) [ppv]\n"
216 "# $7 = tropopause ozone (mean) [ppv]\n"
217 "# $8 = tropopause height (sigma) [km]\n"
218 "# $9 = tropopause pressure (sigma) [hPa]\n"
219 "# $10 = tropopause temperature (sigma) [K]\n"
220 "# $11 = tropopause water vapor (sigma) [ppv]\n"
221 "# $12 = tropopause ozone (sigma) [ppv]\n"
222 "# $13 = number of data points\n"
223 "# $14 = occurrence frequency [%%]\n\n");
224
225 /* Write output... */
226 for (ilat = 0; ilat < nlat; ilat++)
227 fprintf(out, "%.2f %g %g %g %g %g %g %g %g %g %g %g %d %g\n",
228 time0, lats[ilat], zm[ilat], pm[ilat], tm[ilat], qm[ilat],
229 o3m[ilat], zs[ilat], ps[ilat], ts[ilat], qs[ilat], o3s[ilat],
230 n[ilat], 100. * n[ilat] / nt[ilat]);
231
232 /* Close files... */
233 fclose(out);
234
235 return EXIT_SUCCESS;
236}
237
238/*****************************************************************************/
239
241void usage(
242 void) {
243
244 printf("\nMPTRAC tropo_zm tool.\n\n");
245 printf("Extract zonal means from tropopause data.\n");
246 printf("\n");
247 printf("Usage:\n");
248 printf(" tropo_zm <ctl> <zm.tab> <var> <tropo.nc>\n");
249 printf("\n");
250 printf("Arguments:\n");
251 printf(" <ctl> Control file.\n");
252 printf(" <zm.tab> Output table.\n");
253 printf
254 (" <var> Tropopause variable prefix, such as clp, dyn, wmo_1st,\n");
255 printf(" or wmo_2nd.\n");
256 printf(" <tropo.nc> Tropopause netCDF input file.\n");
257 printf("\nFurther information:\n");
258 printf(" Manual: https://slcs-jsc.github.io/mptrac/\n");
259}
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_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
mptrac.h
MPTRAC library declarations.
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: mptrac.h:345
NC
#define NC(cmd)
Execute a NetCDF command and check for errors.
Definition: mptrac.h:1204
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:2102
EY
#define EY
Maximum number of latitudes for meteo data.
Definition: mptrac.h:340
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: mptrac.h:1909
EX
#define EX
Maximum number of longitudes for meteo data.
Definition: mptrac.h:335
SQR
#define SQR(x)
Compute the square of a value.
Definition: mptrac.h:1733
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:2032
NC_GET_DOUBLE
#define NC_GET_DOUBLE(varname, ptr, force)
Retrieve a double-precision variable from a NetCDF file.
Definition: mptrac.h:1263
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:1293
ctl_t
Control parameters.
Definition: mptrac.h:2190
main
int main(int argc, char *argv[])
Definition: tropo_zm.c:46
NT
#define NT
Maximum number of time steps.
Definition: tropo_zm.c:32
usage
void usage(void)
Print command-line help.
Definition: tropo_zm.c:241
Generated by doxygen 1.9.4