51 {
52
54
56
58
60
61 FILE *out;
62
73 zt, tt, plev[
NZ], ps, ts, zs, us, vs, ess, nss, shf, lsm, sst, pbl,
74 pt, pct, pcb, plcl, plfc, pel, cape, cin, o3c, cl, t, u, v, w, pv,
75 h2o, h2ot, o3, lwc, rwc, iwc, swc, cc, lat, lats[
NY], lonm[
NZ][
NY], cw[3];
76
77 static int np[
NZ][
NY], npc[
NZ][
NY], npt[
NZ][
NY], ny, nz, ci[3];
78
79
80 mptrac_alloc(NULL, NULL, &clim, &met, NULL, NULL, NULL, &dd);
81
82
84
85
86 if (argc < 4)
87 ERRMSG(
"Missing or invalid command-line arguments.\n\n"
88 "Usage: met_zm <ctl> <zm.tab> <met0> [<met1> ...]\n\n"
89 "Use -h for full help.");
90
91
93 double z0 =
scan_ctl(argv[1], argc, argv,
"ZM_Z0", -1,
"-999", NULL);
94 double z1 =
scan_ctl(argv[1], argc, argv,
"ZM_Z1", -1,
"-999", NULL);
95 double dz =
scan_ctl(argv[1], argc, argv,
"ZM_DZ", -1,
"-999", NULL);
96 double lon0 =
scan_ctl(argv[1], argc, argv,
"ZM_LON0", -1,
"-360", NULL);
97 double lon1 =
scan_ctl(argv[1], argc, argv,
"ZM_LON1", -1,
"360", NULL);
98 double lat0 =
scan_ctl(argv[1], argc, argv,
"ZM_LAT0", -1,
"-90", NULL);
99 double lat1 =
scan_ctl(argv[1], argc, argv,
"ZM_LAT1", -1,
"90", NULL);
100 double dlat =
scan_ctl(argv[1], argc, argv,
"ZM_DLAT", -1,
"-999", NULL);
101
102
104
105
106 for (int i = 3; i < argc; i++) {
107
108
110 continue;
111
112
113 if (z0 < 0)
115 if (z1 < 0)
116 z1 =
Z(met->
p[met->
np - 1]);
117 nz = 0;
118 if (dz < 0) {
119 for (
int iz = 0; iz < met->
np; iz++)
120 if (
Z(met->
p[iz]) >= z0 &&
Z(met->
p[iz]) <= z1) {
121 plev[nz] = met->
p[iz];
123 ERRMSG(
"Too many pressure levels!");
124 }
125 } else
126 for (z = z0; z <= z1; z += dz) {
129 ERRMSG(
"Too many pressure levels!");
130 }
131
132
133 if (dlat <= 0)
134 dlat = fabs(met->
lat[1] - met->
lat[0]);
135 ny = 0;
136 if (lat0 < -90 && lat1 > 90) {
137 lat0 = gsl_stats_min(met->
lat, 1, (
size_t) met->
ny);
138 lat1 = gsl_stats_max(met->
lat, 1, (
size_t) met->
ny);
139 }
140 for (lat = lat0; lat <= lat1 + 0.001; lat += dlat) {
141 lats[ny] = round(lat * 1e3) / 1e3;
143 ERRMSG(
"Too many latitudes!");
144 }
145
146
147 for (
int ix = 0; ix < met->
nx; ix++)
148 if (met->
lon[ix] >= lon0 && met->
lon[ix] <= lon1)
149 for (int iy = 0; iy < ny; iy++)
150 for (int iz = 0; iz < nz; iz++) {
151
152
154
155
156 timem[iz][iy] += met->
time;
157 lonm[iz][iy] += met->
lon[ix];
158 zm[iz][iy] += z;
159 tm[iz][iy] += t;
160 um[iz][iy] += u;
161 vm[iz][iy] += v;
162 wm[iz][iy] += w;
163 pvm[iz][iy] += pv;
164 h2om[iz][iy] += h2o;
165 o3m[iz][iy] += o3;
166 lwcm[iz][iy] += lwc;
167 rwcm[iz][iy] += rwc;
168 iwcm[iz][iy] += iwc;
169 swcm[iz][iy] += swc;
170 ccm[iz][iy] += cc;
171 psm[iz][iy] += ps;
172 tsm[iz][iy] += ts;
173 zsm[iz][iy] += zs;
174 usm[iz][iy] += us;
175 vsm[iz][iy] += vs;
176 essm[iz][iy] += ess;
177 nssm[iz][iy] += nss;
178 shfm[iz][iy] += shf;
179 lsmm[iz][iy] += lsm;
180 sstm[iz][iy] += sst;
181 pblm[iz][iy] += pbl;
182 pctm[iz][iy] += pct;
183 pcbm[iz][iy] += pcb;
184 clm[iz][iy] += cl;
185 if (isfinite(plfc) && isfinite(pel) && cape >= ctl.
conv_cape
187 plclm[iz][iy] += plcl;
188 plfcm[iz][iy] += plfc;
189 pelm[iz][iy] += pel;
190 capem[iz][iy] += cape;
191 cinm[iz][iy] += cin;
192 npc[iz][iy]++;
193 }
194 if (isfinite(pt)) {
195 ptm[iz][iy] += pt;
196 ztm[iz][iy] += zt;
197 ttm[iz][iy] += tt;
198 h2otm[iz][iy] += h2ot;
199 npt[iz][iy]++;
200 }
201 o3cm[iz][iy] += o3c;
202 rhm[iz][iy] +=
RH(plev[iz], t, h2o);
203 rhicem[iz][iy] +=
RHICE(plev[iz], t, h2o);
204 tdewm[iz][iy] +=
TDEW(plev[iz], h2o);
205 ticem[iz][iy] +=
TICE(plev[iz], h2o);
206 hno3m[iz][iy] +=
208 tnatm[iz][iy] +=
211 plev[iz]));
212 ohm[iz][iy] +=
214 plev[iz]);
215 h2o2m[iz][iy]
217 ho2m[iz][iy]
219 o1dm[iz][iy]
221 np[iz][iy]++;
222 }
223 }
224
225
226 LOG(1,
"Write meteorological data file: %s", argv[2]);
227 if (!(out = fopen(argv[2], "w")))
228 ERRMSG(
"Cannot create file!");
229
230
232
233
234 for (int iz = 0; iz < nz; iz++) {
235 fprintf(out, "\n");
236 for (int iy = 0; iy < ny; iy++)
237 fprintf(out,
238 "%.2f %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
239 " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
240 " %g %g %g %g %g %g %g %g %g %g %g %g %g %d %d %d\n",
241 timem[iz][iy] / np[iz][iy],
Z(plev[iz]),
242 lonm[iz][iy] / np[iz][iy], lats[iy],
243 plev[iz], tm[iz][iy] / np[iz][iy], um[iz][iy] / np[iz][iy],
244 vm[iz][iy] / np[iz][iy], wm[iz][iy] / np[iz][iy],
245 h2om[iz][iy] / np[iz][iy], o3m[iz][iy] / np[iz][iy],
246 zm[iz][iy] / np[iz][iy], pvm[iz][iy] / np[iz][iy],
247 psm[iz][iy] / np[iz][iy], tsm[iz][iy] / np[iz][iy],
248 zsm[iz][iy] / np[iz][iy], usm[iz][iy] / np[iz][iy],
249 vsm[iz][iy] / np[iz][iy], essm[iz][iy] / np[iz][iy],
250 nssm[iz][iy] / np[iz][iy], shfm[iz][iy] / np[iz][iy],
251 lsmm[iz][iy] / np[iz][iy],
252 sstm[iz][iy] / np[iz][iy], ptm[iz][iy] / npt[iz][iy],
253 ztm[iz][iy] / npt[iz][iy], ttm[iz][iy] / npt[iz][iy],
254 h2otm[iz][iy] / npt[iz][iy], lwcm[iz][iy] / np[iz][iy],
255 rwcm[iz][iy] / np[iz][iy], iwcm[iz][iy] / np[iz][iy],
256 swcm[iz][iy] / np[iz][iy], ccm[iz][iy] / np[iz][iy],
257 clm[iz][iy] / np[iz][iy], pctm[iz][iy] / np[iz][iy],
258 pcbm[iz][iy] / np[iz][iy], plclm[iz][iy] / npc[iz][iy],
259 plfcm[iz][iy] / npc[iz][iy], pelm[iz][iy] / npc[iz][iy],
260 capem[iz][iy] / npc[iz][iy], cinm[iz][iy] / npc[iz][iy],
261 rhm[iz][iy] / np[iz][iy], rhicem[iz][iy] / np[iz][iy],
262 tdewm[iz][iy] / np[iz][iy], ticem[iz][iy] / np[iz][iy],
263 tnatm[iz][iy] / np[iz][iy], hno3m[iz][iy] / np[iz][iy],
264 ohm[iz][iy] / np[iz][iy], h2o2m[iz][iy] / np[iz][iy],
265 ho2m[iz][iy] / np[iz][iy], o1dm[iz][iy] / np[iz][iy],
266 pblm[iz][iy] / np[iz][iy], o3cm[iz][iy] / np[iz][iy],
267 np[iz][iy], npt[iz][iy], npc[iz][iy]);
268 }
269
270
271 fclose(out);
272
273
274 mptrac_free(NULL, NULL, clim, met, NULL, NULL, NULL, dd);
275
276 return EXIT_SUCCESS;
277}
#define NY
Maximum number of latitudes.
#define NZ
Maximum number of altitudes.
double clim_zm(const clim_zm_t *zm, const double t, const double lat, const double p)
Interpolates monthly mean zonal mean climatological variables.
void mptrac_free(ctl_t *ctl, cache_t *cache, clim_t *clim, met_t *met0, met_t *met1, atm_t *atm, depo_t *depo, dd_t *dd)
Frees memory resources allocated for MPTRAC.
double nat_temperature(const double p, const double h2o, const double hno3)
Calculates the nitric acid trihydrate (NAT) temperature.
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.
void mptrac_read_clim(const ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
double clim_oh(const ctl_t *ctl, const clim_t *clim, const double t, const double lon, const double lat, const double p)
Calculates the hydroxyl radical (OH) concentration from climatology data, with an optional diurnal co...
void mptrac_alloc(ctl_t **ctl, cache_t **cache, clim_t **clim, met_t **met0, met_t **met1, atm_t **atm, depo_t **depo, dd_t **dd)
Allocates and initializes memory resources for MPTRAC.
int mptrac_read_met(const char *filename, const ctl_t *ctl, const clim_t *clim, met_t *met, dd_t *dd)
Reads meteorological data from a file, supporting multiple formats and MPI broadcasting.
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.
#define INTPOL_SPACE_ALL(p, lon, lat)
Interpolate multiple meteorological variables in space.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
#define USAGE
Print usage information on -h or --help.
#define Z(p)
Convert pressure to altitude.
#define P(z)
Compute pressure at given altitude.
#define MET_HEADER
Write header for meteorological data file.
#define TICE(p, h2o)
Calculate frost point temperature (WMO, 2018).
#define RHICE(p, t, h2o)
Compute relative humidity over ice.
#define RH(p, t, h2o)
Compute relative humidity over water.
#define LOG(level,...)
Print a log message with a specified logging level.
#define TDEW(p, h2o)
Calculate dew point temperature.
clim_zm_t ho2
HO2 zonal means.
clim_zm_t hno3
HNO3 zonal means.
clim_zm_t o1d
O(1D) zonal means.
clim_zm_t h2o2
H2O2 zonal means.
double conv_cape
CAPE threshold for convection module [J/kg].
double conv_cin
CIN threshold for convection module [J/kg].
Domain decomposition data structure.
int nx
Number of longitudes.
int ny
Number of latitudes.
int np
Number of pressure levels.
double lon[EX]
Longitudes [deg].
double lat[EY]
Latitudes [deg].
double p[EP]
Pressure levels [hPa].