MPTRAC
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-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 atm_t *atm;
45
46 static FILE *out;
47
48 static char varname[LEN];
49
50 static double times[NT], lons[EX], lats[EY], time0, time1, z0, z0sig, p0,
51 p0sig, t0, t0sig, q0, q0sig, o30, o30sig;
52
53 static float help[EX * EY], tropo_z0[EX][EY], tropo_z1[EX][EY],
54 tropo_p0[EX][EY], tropo_p1[EX][EY], tropo_t0[EX][EY], tropo_t1[EX][EY],
55 tropo_q0[EX][EY], tropo_q1[EX][EY], tropo_o30[EX][EY], tropo_o31[EX][EY];
56
57 static int it_old =
58 -999, ncid, varid, varid_z, varid_p, varid_t, varid_q, varid_o3, h2o, o3,
59 ntime, nlon, nlat, ilon, ilat;
60
61 static size_t count[10], start[10];
62
63 /* Allocate... */
64 ALLOC(atm, atm_t, 1);
65
66 /* Check arguments... */
67 if (argc < 5)
68 ERRMSG("Give parameters: <ctl> <sample.tab> <tropo.nc> <var> <atm_in>");
69
70 /* Read control parameters... */
71 read_ctl(argv[1], argc, argv, &ctl);
72 int method =
73 (int) scan_ctl(argv[1], argc, argv, "TROPO_SAMPLE_METHOD", -1, "1", NULL);
74
75 /* Read atmospheric data... */
76 if (!read_atm(argv[5], &ctl, atm))
77 ERRMSG("Cannot open file!");
78
79 /* Open tropopause file... */
80 LOG(1, "Read tropopause data: %s", argv[3]);
81 if (nc_open(argv[3], NC_NOWRITE, &ncid) != NC_NOERR)
82 ERRMSG("Cannot open file!");
83
84 /* Get dimensions... */
85 NC_INQ_DIM("time", &ntime, 1, NT);
86 NC_INQ_DIM("lat", &nlat, 1, EY);
87 NC_INQ_DIM("lon", &nlon, 1, EX);
88
89 /* Read coordinates... */
90 NC_GET_DOUBLE("time", times, 1);
91 NC_GET_DOUBLE("lat", lats, 1);
92 NC_GET_DOUBLE("lon", lons, 1);
93
94 /* Get variable indices... */
95 sprintf(varname, "%s_z", argv[4]);
96 NC(nc_inq_varid(ncid, varname, &varid_z));
97 sprintf(varname, "%s_p", argv[4]);
98 NC(nc_inq_varid(ncid, varname, &varid_p));
99 sprintf(varname, "%s_t", argv[4]);
100 NC(nc_inq_varid(ncid, varname, &varid_t));
101 sprintf(varname, "%s_q", argv[4]);
102 h2o = (nc_inq_varid(ncid, varname, &varid_q) == NC_NOERR);
103 sprintf(varname, "%s_o3", argv[4]);
104 o3 = (nc_inq_varid(ncid, varname, &varid_o3) == NC_NOERR);
105
106 /* Set dimensions... */
107 count[0] = 1;
108 count[1] = (size_t) nlat;
109 count[2] = (size_t) nlon;
110
111 /* Create file... */
112 LOG(1, "Write tropopause sample data: %s", argv[2]);
113 if (!(out = fopen(argv[2], "w")))
114 ERRMSG("Cannot create file!");
115
116 /* Write header... */
117 fprintf(out,
118 "# $1 = time [s]\n"
119 "# $2 = altitude [km]\n"
120 "# $3 = longitude [deg]\n" "# $4 = latitude [deg]\n");
121 for (int iq = 0; iq < ctl.nq; iq++)
122 fprintf(out, "# $%i = %s [%s]\n", iq + 5, ctl.qnt_name[iq],
123 ctl.qnt_unit[iq]);
124 fprintf(out, "# $%d = tropopause height (mean) [km]\n", 5 + ctl.nq);
125 fprintf(out, "# $%d = tropopause pressure (mean) [hPa]\n", 6 + ctl.nq);
126 fprintf(out, "# $%d = tropopause temperature (mean) [K]\n", 7 + ctl.nq);
127 fprintf(out, "# $%d = tropopause water vapor (mean) [ppv]\n", 8 + ctl.nq);
128 fprintf(out, "# $%d = tropopause ozone (mean) [ppv]\n", 9 + ctl.nq);
129 fprintf(out, "# $%d = tropopause height (sigma) [km]\n", 10 + ctl.nq);
130 fprintf(out, "# $%d = tropopause pressure (sigma) [hPa]\n", 11 + ctl.nq);
131 fprintf(out, "# $%d = tropopause temperature (sigma) [K]\n", 12 + ctl.nq);
132 fprintf(out, "# $%d = tropopause water vapor (sigma) [ppv]\n", 13 + ctl.nq);
133 fprintf(out, "# $%d = tropopause ozone (sigma) [ppv]\n\n", 14 + ctl.nq);
134
135 /* Loop over particles... */
136 for (int ip = 0; ip < atm->np; ip++) {
137
138 /* Check temporal ordering... */
139 if (ip > 0 && atm->time[ip] < atm->time[ip - 1])
140 ERRMSG("Time must be ascending!");
141
142 /* Check range... */
143 if (atm->time[ip] < times[0] || atm->time[ip] > times[ntime - 1])
144 continue;
145
146 /* Read data... */
147 int it = locate_irr(times, (int) ntime, atm->time[ip]);
148 if (it != it_old) {
149
150 time0 = times[it];
151 start[0] = (size_t) it;
152 NC(nc_get_vara_float(ncid, varid_z, start, count, help));
153 for (ilon = 0; ilon < nlon; ilon++)
154 for (ilat = 0; ilat < nlat; ilat++)
155 tropo_z0[ilon][ilat] = help[ilat * nlon + ilon];
156 NC(nc_get_vara_float(ncid, varid_p, start, count, help));
157 for (ilon = 0; ilon < nlon; ilon++)
158 for (ilat = 0; ilat < nlat; ilat++)
159 tropo_p0[ilon][ilat] = help[ilat * nlon + ilon];
160 NC(nc_get_vara_float(ncid, varid_t, start, count, help));
161 for (ilon = 0; ilon < nlon; ilon++)
162 for (ilat = 0; ilat < nlat; ilat++)
163 tropo_t0[ilon][ilat] = help[ilat * nlon + ilon];
164 if (h2o) {
165 NC(nc_get_vara_float(ncid, varid_q, start, count, help));
166 for (ilon = 0; ilon < nlon; ilon++)
167 for (ilat = 0; ilat < nlat; ilat++)
168 tropo_q0[ilon][ilat] = help[ilat * nlon + ilon];
169 } else
170 for (ilon = 0; ilon < nlon; ilon++)
171 for (ilat = 0; ilat < nlat; ilat++)
172 tropo_q0[ilon][ilat] = NAN;
173 if (o3) {
174 NC(nc_get_vara_float(ncid, varid_o3, start, count, help));
175 for (ilon = 0; ilon < nlon; ilon++)
176 for (ilat = 0; ilat < nlat; ilat++)
177 tropo_o30[ilon][ilat] = help[ilat * nlon + ilon];
178 } else
179 for (ilon = 0; ilon < nlon; ilon++)
180 for (ilat = 0; ilat < nlat; ilat++)
181 tropo_o30[ilon][ilat] = NAN;
182
183 time1 = times[it + 1];
184 start[0] = (size_t) it + 1;
185 NC(nc_get_vara_float(ncid, varid_z, start, count, help));
186 for (ilon = 0; ilon < nlon; ilon++)
187 for (ilat = 0; ilat < nlat; ilat++)
188 tropo_z1[ilon][ilat] = help[ilat * nlon + ilon];
189 NC(nc_get_vara_float(ncid, varid_p, start, count, help));
190 for (ilon = 0; ilon < nlon; ilon++)
191 for (ilat = 0; ilat < nlat; ilat++)
192 tropo_p1[ilon][ilat] = help[ilat * nlon + ilon];
193 NC(nc_get_vara_float(ncid, varid_t, start, count, help));
194 for (ilon = 0; ilon < nlon; ilon++)
195 for (ilat = 0; ilat < nlat; ilat++)
196 tropo_t1[ilon][ilat] = help[ilat * nlon + ilon];
197 if (h2o) {
198 NC(nc_get_vara_float(ncid, varid_q, start, count, help));
199 for (ilon = 0; ilon < nlon; ilon++)
200 for (ilat = 0; ilat < nlat; ilat++)
201 tropo_q1[ilon][ilat] = help[ilat * nlon + ilon];
202 } else
203 for (ilon = 0; ilon < nlon; ilon++)
204 for (ilat = 0; ilat < nlat; ilat++)
205 tropo_q1[ilon][ilat] = NAN;
206 if (o3) {
207 NC(nc_get_vara_float(ncid, varid_o3, start, count, help));
208 for (ilon = 0; ilon < nlon; ilon++)
209 for (ilat = 0; ilat < nlat; ilat++)
210 tropo_o31[ilon][ilat] = help[ilat * nlon + ilon];
211 } else
212 for (ilon = 0; ilon < nlon; ilon++)
213 for (ilat = 0; ilat < nlat; ilat++)
214 tropo_o31[ilon][ilat] = NAN;
215 }
216 it_old = it;
217
218 /* Interpolate... */
219 intpol_tropo_3d(time0, tropo_z0, time1, tropo_z1,
220 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
221 atm->lat[ip], method, &z0, &z0sig);
222 intpol_tropo_3d(time0, tropo_p0, time1, tropo_p1,
223 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
224 atm->lat[ip], method, &p0, &p0sig);
225 intpol_tropo_3d(time0, tropo_t0, time1, tropo_t1,
226 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
227 atm->lat[ip], method, &t0, &t0sig);
228 intpol_tropo_3d(time0, tropo_q0, time1, tropo_q1,
229 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
230 atm->lat[ip], method, &q0, &q0sig);
231 intpol_tropo_3d(time0, tropo_o30, time1, tropo_o31,
232 lons, lats, nlon, nlat, atm->time[ip], atm->lon[ip],
233 atm->lat[ip], method, &o30, &o30sig);
234
235 /* Write output... */
236 fprintf(out, "%.2f %g %g %g", atm->time[ip], Z(atm->p[ip]),
237 atm->lon[ip], atm->lat[ip]);
238 for (int iq = 0; iq < ctl.nq; iq++) {
239 fprintf(out, " ");
240 fprintf(out, ctl.qnt_format[iq], atm->q[iq][ip]);
241 }
242 fprintf(out, " %g %g %g %g %g %g %g %g %g %g\n",
243 z0, p0, t0, q0, o30, z0sig, p0sig, t0sig, q0sig, o30sig);
244 }
245
246 /* Close files... */
247 fclose(out);
248 NC(nc_close(ncid));
249
250 /* Free... */
251 free(atm);
252
253 return EXIT_SUCCESS;
254}
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:2116
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:8098
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:1738
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
read_atm
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.
Definition: mptrac.c:4215
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
Z
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:1726
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
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
atm_t
Air parcel data.
Definition: mptrac.h:3120
atm_t::time
double time[NP]
Time [s].
Definition: mptrac.h:3126
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3135
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3132
atm_t::np
int np
Number of air parcels.
Definition: mptrac.h:3123
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3138
atm_t::p
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3129
ctl_t
Control parameters.
Definition: mptrac.h:2155
ctl_t::qnt_format
char qnt_format[NQ][LEN]
Quantity output format.
Definition: mptrac.h:2192
ctl_t::qnt_unit
char qnt_unit[NQ][LEN]
Quantity units.
Definition: mptrac.h:2189
ctl_t::qnt_name
char qnt_name[NQ][LEN]
Quantity names.
Definition: mptrac.h:2183
ctl_t::nq
int nq
Number of quantities.
Definition: mptrac.h:2180
main
int main(int argc, char *argv[])
Definition: tropo_sample.c:38
NT
#define NT
Maximum number of time steps.
Definition: tropo_sample.c:32
Generated on Thu Nov 21 2024 13:53:30 for MPTRAC by doxygen 1.9.4