MPTRAC
tropo_clim.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-2024 Forschungszentrum Juelich GmbH
18*/
19
25#include "mptrac.h"
26
27/* ------------------------------------------------------------
28 Dimensions...
29 ------------------------------------------------------------ */
30
32#define NT 744
33
34/* ------------------------------------------------------------
35 Main...
36 ------------------------------------------------------------ */
37
38int main(
39 int argc,
40 char *argv[]) {
41
42 ctl_t ctl;
43
44 static FILE *out;
45
46 static char varname[LEN];
47
48 static double lons[EX], lats[EY], *zm, *zs, *pm, *ps, *tm, *ts, *qm, *qs,
49 *o3m, *o3s;
50
51 static float *tropo_z0, *tropo_p0, *tropo_t0, *tropo_q0, *tropo_o30;
52
53 static int ncid, varid, varid_z, varid_p, varid_t, varid_q, varid_o3, h2o,
54 o3, *n, *nt, ntime, nlon, nlat;
55
56 static size_t count[10], start[10];
57
58 /* Allocate... */
59 ALLOC(zm, double,
60 EX * EY);
61 ALLOC(zs, double,
62 EX * EY);
63 ALLOC(pm, double,
64 EX * EY);
65 ALLOC(ps, double,
66 EX * EY);
67 ALLOC(tm, double,
68 EX * EY);
69 ALLOC(ts, double,
70 EX * EY);
71 ALLOC(qm, double,
72 EX * EY);
73 ALLOC(qs, double,
74 EX * EY);
75 ALLOC(o3m, double,
76 EX * EY);
77 ALLOC(o3s, double,
78 EX * EY);
79
80 ALLOC(tropo_z0, float,
81 EX * EY);
82 ALLOC(tropo_p0, float,
83 EX * EY);
84 ALLOC(tropo_t0, float,
85 EX * EY);
86 ALLOC(tropo_q0, float,
87 EX * EY);
88 ALLOC(tropo_o30, float,
89 EX * EY);
90
91 ALLOC(n, int,
92 EX * EY);
93 ALLOC(nt, int,
94 EX * EY);
95
96 /* Check arguments... */
97 if (argc < 5)
98 ERRMSG("Give parameters: <ctl> <clim.tab> <var> <tropo.nc>");
99
100 /* Read control parameters... */
101 read_ctl(argv[1], argc, argv, &ctl);
102
103 /* Loop over tropopause files... */
104 for (int iarg = 4; iarg < argc; iarg++) {
105
106 /* Open tropopause file... */
107 LOG(1, "Read tropopause data: %s", argv[iarg]);
108 if (nc_open(argv[iarg], NC_NOWRITE, &ncid) != NC_NOERR)
109 ERRMSG("Cannot open file!");
110
111 /* Get dimensions... */
112 NC_INQ_DIM("time", &ntime, 1, NT);
113 NC_INQ_DIM("lat", &nlat, 1, EY);
114 NC_INQ_DIM("lon", &nlon, 1, EX);
115
116 /* Read coordinates... */
117 NC_GET_DOUBLE("lat", lats, 1);
118 NC_GET_DOUBLE("lon", lons, 1);
119
120 /* Get variable indices... */
121 sprintf(varname, "%s_z", argv[3]);
122 NC(nc_inq_varid(ncid, varname, &varid_z));
123 sprintf(varname, "%s_p", argv[3]);
124 NC(nc_inq_varid(ncid, varname, &varid_p));
125 sprintf(varname, "%s_t", argv[3]);
126 NC(nc_inq_varid(ncid, varname, &varid_t));
127 sprintf(varname, "%s_q", argv[3]);
128 h2o = (nc_inq_varid(ncid, varname, &varid_q) == NC_NOERR);
129 sprintf(varname, "%s_o3", argv[3]);
130 o3 = (nc_inq_varid(ncid, varname, &varid_o3) == NC_NOERR);
131
132 /* Set dimensions... */
133 count[0] = 1;
134 count[1] = (size_t) nlat;
135 count[2] = (size_t) nlon;
136
137 /* Loop over time steps... */
138 for (int it = 0; it < ntime; it++) {
139
140 /* Read data... */
141 start[0] = (size_t) it;
142 NC(nc_get_vara_float(ncid, varid_z, start, count, tropo_z0));
143 NC(nc_get_vara_float(ncid, varid_p, start, count, tropo_p0));
144 NC(nc_get_vara_float(ncid, varid_t, start, count, tropo_t0));
145 if (h2o) {
146 NC(nc_get_vara_float(ncid, varid_q, start, count, tropo_q0));
147 } else
148 for (int i = 0; i < nlon * nlat; i++)
149 tropo_q0[i] = NAN;
150 if (o3) {
151 NC(nc_get_vara_float(ncid, varid_o3, start, count, tropo_o30));
152 } else
153 for (int i = 0; i < nlon * nlat; i++)
154 tropo_o30[i] = NAN;
155
156 /* Averaging... */
157 for (int i = 0; i < nlon * nlat; i++) {
158 nt[i]++;
159 if (isfinite(tropo_z0[i])
160 && isfinite(tropo_p0[i])
161 && isfinite(tropo_t0[i])
162 && (!h2o || isfinite(tropo_q0[i]))
163 && (!o3 || isfinite(tropo_o30[i]))) {
164 zm[i] += tropo_z0[i];
165 zs[i] += SQR(tropo_z0[i]);
166 pm[i] += tropo_p0[i];
167 ps[i] += SQR(tropo_p0[i]);
168 tm[i] += tropo_t0[i];
169 ts[i] += SQR(tropo_t0[i]);
170 qm[i] += tropo_q0[i];
171 qs[i] += SQR(tropo_q0[i]);
172 o3m[i] += tropo_o30[i];
173 o3s[i] += SQR(tropo_o30[i]);
174 n[i]++;
175 }
176 }
177 }
178
179 /* Close files... */
180 NC(nc_close(ncid));
181 }
182
183 /* Normalize... */
184 for (int i = 0; i < nlon * nlat; i++)
185 if (n[i] > 0) {
186 zm[i] /= n[i];
187 pm[i] /= n[i];
188 tm[i] /= n[i];
189 qm[i] /= n[i];
190 o3m[i] /= n[i];
191 double aux = zs[i] / n[i] - SQR(zm[i]);
192 zs[i] = aux > 0 ? sqrt(aux) : 0.0;
193 aux = ps[i] / n[i] - SQR(pm[i]);
194 ps[i] = aux > 0 ? sqrt(aux) : 0.0;
195 aux = ts[i] / n[i] - SQR(tm[i]);
196 ts[i] = aux > 0 ? sqrt(aux) : 0.0;
197 aux = qs[i] / n[i] - SQR(qm[i]);
198 qs[i] = aux > 0 ? sqrt(aux) : 0.0;
199 aux = o3s[i] / n[i] - SQR(o3m[i]);
200 o3s[i] = aux > 0 ? sqrt(aux) : 0.0;
201 }
202
203 /* Create file... */
204 LOG(1, "Write tropopause climatological 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 = longitude [deg]\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");
224
225 /* Write output... */
226 for (int ilat = 0; ilat < nlat; ilat++) {
227 fprintf(out, "\n");
228 for (int ilon = 0; ilon < nlon; ilon++)
229 fprintf(out, "%g %g %g %g %g %g %g %g %g %g %g %g %d %g\n",
230 lons[ilon], lats[ilat], zm[ilat * nlon + ilon],
231 pm[ilat * nlon + ilon], tm[ilat * nlon + ilon],
232 qm[ilat * nlon + ilon], o3m[ilat * nlon + ilon],
233 zs[ilat * nlon + ilon], ps[ilat * nlon + ilon],
234 ts[ilat * nlon + ilon], qs[ilat * nlon + ilon],
235 o3s[ilat * nlon + ilon], n[ilat * nlon + ilon],
236 100. * n[ilat * nlon + ilon] / nt[ilat * nlon + ilon]);
237 }
238
239 /* Close files... */
240 fclose(out);
241
242 /* Free... */
243 free(zm);
244 free(zs);
245 free(pm);
246 free(ps);
247 free(tm);
248 free(ts);
249 free(qm);
250 free(qs);
251 free(o3m);
252 free(o3s);
253
254 free(tropo_z0);
255 free(tropo_p0);
256 free(tropo_t0);
257 free(tropo_q0);
258 free(tropo_o30);
259
260 free(n);
261 free(nt);
262
263 return EXIT_SUCCESS;
264}
read_ctl
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:4805
mptrac.h
MPTRAC library declarations.
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: mptrac.h:236
NC
#define NC(cmd)
Execute a NetCDF command and check for errors.
Definition: mptrac.h:1001
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:1901
EY
#define EY
Maximum number of latitudes for meteo data.
Definition: mptrac.h:269
EX
#define EX
Maximum number of longitudes for meteo data.
Definition: mptrac.h:262
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:347
SQR
#define SQR(x)
Compute the square of a value.
Definition: mptrac.h:1542
NC_INQ_DIM
#define NC_INQ_DIM(dimname, ptr, min, max)
Inquire the length of a dimension in a NetCDF file.
Definition: mptrac.h:1088
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:1831
NC_GET_DOUBLE
#define NC_GET_DOUBLE(varname, ptr, force)
Retrieve a double-precision variable from a NetCDF file.
Definition: mptrac.h:1060
ctl_t
Control parameters.
Definition: mptrac.h:2155
main
int main(int argc, char *argv[])
Definition: tropo_clim.c:38
NT
#define NT
Maximum number of time steps.
Definition: tropo_clim.c:32
Generated on Thu Nov 21 2024 13:53:30 for MPTRAC by doxygen 1.9.4