51 {
52
54
56
58
60
61 FILE *out;
62
63 static double *timem, ps, *psm, ts, *tsm, zs, *zsm, us, *usm, vs, *vsm,
64 ess, *essm, nss, *nssm, shf, *shfm, lsm, *lsmm, sst, *sstm, pbl, *pblm,
65 pt, *ptm, t, *pm, *tm, u, *um, v, *vm, w, *wm, h2o, *h2om, h2ot, *h2otm,
66 o3, *o3m, *hno3m, *ohm, *h2o2m, *ho2m, *o1dm, *tdewm, *ticem, *tnatm,
67 lwc, *lwcm, rwc, *rwcm, iwc, *iwcm, swc, *swcm, cc, *ccm, z, *zm,
68 pv, *pvm, zt, *ztm, tt, *ttm, pct, *pctm, pcb, *pcbm, cl, *clm,
69 plcl, *plclm, plfc, *plfcm, pel, *pelm, cape, *capem, cin, *cinm,
70 o3c, *o3cm, *rhm, *rhicem, ptop, pbot, t0, lon, lons[
NX], lat, lats[
NY];
71
72 static int *np, *npc, *npt, nx, ny;
73
74
75 mptrac_alloc(NULL, NULL, &clim, &met, NULL, NULL, NULL, &dd);
172 ALLOC(rhicem,
double,
180
181
183
184
185 if (argc < 4)
186 ERRMSG(
"Missing or invalid command-line arguments.\n\n"
187 "Usage: met_map <ctl> <map.tab> <met0> [<met1> ...]\n\n"
188 "Use -h for full help.");
189
190
192 double p0 =
P(
scan_ctl(argv[1], argc, argv,
"MAP_Z0", -1,
"10", NULL));
193 double lon0 =
scan_ctl(argv[1], argc, argv,
"MAP_LON0", -1,
"-180", NULL);
194 double lon1 =
scan_ctl(argv[1], argc, argv,
"MAP_LON1", -1,
"180", NULL);
195 double dlon =
scan_ctl(argv[1], argc, argv,
"MAP_DLON", -1,
"-999", NULL);
196 double lat0 =
scan_ctl(argv[1], argc, argv,
"MAP_LAT0", -1,
"-90", NULL);
197 double lat1 =
scan_ctl(argv[1], argc, argv,
"MAP_LAT1", -1,
"90", NULL);
198 double dlat =
scan_ctl(argv[1], argc, argv,
"MAP_DLAT", -1,
"-999", NULL);
199 double theta =
scan_ctl(argv[1], argc, argv,
"MAP_THETA", -1,
"-999", NULL);
200
201
203
204
205 for (int i = 3; i < argc; i++) {
206
207
209 continue;
210
211
212 if (dlon <= 0)
213 dlon = fabs(met->
lon[1] - met->
lon[0]);
214 if (dlat <= 0)
215 dlat = fabs(met->
lat[1] - met->
lat[0]);
216 if (lon0 < -360 && lon1 > 360) {
217 lon0 = gsl_stats_min(met->
lon, 1, (
size_t) met->
nx);
218 lon1 = gsl_stats_max(met->
lon, 1, (
size_t) met->
nx);
219 }
220 nx = ny = 0;
221 for (lon = lon0; lon <= lon1 + 0.001; lon += dlon) {
222 lons[nx] = round(lon * 1e3) / 1e3;
224 ERRMSG(
"Too many longitudes!");
225 }
226 if (lat0 < -90 && lat1 > 90) {
227 lat0 = gsl_stats_min(met->
lat, 1, (
size_t) met->
ny);
228 lat1 = gsl_stats_max(met->
lat, 1, (
size_t) met->
ny);
229 }
230 for (lat = lat0; lat <= lat1 + 0.001; lat += dlat) {
231 lats[ny] = round(lat * 1e3) / 1e3;
233 ERRMSG(
"Too many latitudes!");
234 }
235
236
237 for (int ix = 0; ix < nx; ix++)
238 for (int iy = 0; iy < ny; iy++) {
239
240
242 if (theta > 0) {
243 ptop = met->
p[met->
np - 1];
245 do {
246 p0 = 0.5 * (ptop + pbot);
248 &t0, ci, cw, 1);
249 if (
THETA(p0, t0) > theta)
250 ptop = p0;
251 else
252 pbot = p0;
253 } while (fabs(ptop - pbot) > 1e-5);
254 }
255
256
258
259
260 timem[iy * nx + ix] += met->
time;
261 zm[iy * nx + ix] += z;
262 pm[iy * nx + ix] += p0;
263 tm[iy * nx + ix] += t;
264 um[iy * nx + ix] += u;
265 vm[iy * nx + ix] += v;
266 wm[iy * nx + ix] += w;
267 pvm[iy * nx + ix] += pv;
268 h2om[iy * nx + ix] += h2o;
269 o3m[iy * nx + ix] += o3;
270 lwcm[iy * nx + ix] += lwc;
271 rwcm[iy * nx + ix] += rwc;
272 iwcm[iy * nx + ix] += iwc;
273 swcm[iy * nx + ix] += swc;
274 ccm[iy * nx + ix] += cc;
275 psm[iy * nx + ix] += ps;
276 tsm[iy * nx + ix] += ts;
277 zsm[iy * nx + ix] += zs;
278 usm[iy * nx + ix] += us;
279 vsm[iy * nx + ix] += vs;
280 essm[iy * nx + ix] += ess;
281 nssm[iy * nx + ix] += nss;
282 shfm[iy * nx + ix] += shf;
283 lsmm[iy * nx + ix] += lsm;
284 sstm[iy * nx + ix] += sst;
285 pblm[iy * nx + ix] += pbl;
286 pctm[iy * nx + ix] += pct;
287 pcbm[iy * nx + ix] += pcb;
288 clm[iy * nx + ix] += cl;
289 if (isfinite(plfc) && isfinite(pel) && cape >= ctl.
conv_cape
291 plclm[iy * nx + ix] += plcl;
292 plfcm[iy * nx + ix] += plfc;
293 pelm[iy * nx + ix] += pel;
294 capem[iy * nx + ix] += cape;
295 cinm[iy * nx + ix] += cin;
296 npc[iy * nx + ix]++;
297 }
298 if (isfinite(pt)) {
299 ptm[iy * nx + ix] += pt;
300 ztm[iy * nx + ix] += zt;
301 ttm[iy * nx + ix] += tt;
302 h2otm[iy * nx + ix] += h2ot;
303 npt[iy * nx + ix]++;
304 }
305 o3cm[iy * nx + ix] += o3c;
307 tnatm[iy * nx + ix] +=
310 ohm[iy * nx + ix] +=
311 clim_oh(&ctl, clim, met->
time, lons[ix], lats[iy], p0);
313 ho2m[iy * nx + ix] +=
clim_zm(&clim->
ho2, met->
time, lats[iy], p0);
314 o1dm[iy * nx + ix] +=
clim_zm(&clim->
o1d, met->
time, lats[iy], p0);
315 rhm[iy * nx + ix] +=
RH(p0, t, h2o);
316 rhicem[iy * nx + ix] +=
RHICE(p0, t, h2o);
317 tdewm[iy * nx + ix] +=
TDEW(p0, h2o);
318 ticem[iy * nx + ix] +=
TICE(p0, h2o);
319 np[iy * nx + ix]++;
320 }
321 }
322
323
324 LOG(1,
"Write meteorological data file: %s", argv[2]);
325 if (!(out = fopen(argv[2], "w")))
326 ERRMSG(
"Cannot create file!");
327
328
330
331
332 for (int iy = 0; iy < ny; iy++) {
333 fprintf(out, "\n");
334 for (int ix = 0; ix < nx; ix++)
335 fprintf(out,
336 "%.2f %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
337 " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
338 " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %d %d %d\n",
339 timem[iy * nx + ix] / np[iy * nx + ix],
340 Z(pm[iy * nx + ix] / np[iy * nx + ix]), lons[ix], lats[iy],
341 pm[iy * nx + ix] / np[iy * nx + ix],
342 tm[iy * nx + ix] / np[iy * nx + ix],
343 um[iy * nx + ix] / np[iy * nx + ix],
344 vm[iy * nx + ix] / np[iy * nx + ix],
345 wm[iy * nx + ix] / np[iy * nx + ix],
346 h2om[iy * nx + ix] / np[iy * nx + ix],
347 o3m[iy * nx + ix] / np[iy * nx + ix],
348 zm[iy * nx + ix] / np[iy * nx + ix],
349 pvm[iy * nx + ix] / np[iy * nx + ix],
350 psm[iy * nx + ix] / np[iy * nx + ix],
351 tsm[iy * nx + ix] / np[iy * nx + ix],
352 zsm[iy * nx + ix] / np[iy * nx + ix],
353 usm[iy * nx + ix] / np[iy * nx + ix],
354 vsm[iy * nx + ix] / np[iy * nx + ix],
355 essm[iy * nx + ix] / np[iy * nx + ix],
356 nssm[iy * nx + ix] / np[iy * nx + ix],
357 shfm[iy * nx + ix] / np[iy * nx + ix],
358 lsmm[iy * nx + ix] / np[iy * nx + ix],
359 sstm[iy * nx + ix] / np[iy * nx + ix],
360 ptm[iy * nx + ix] / npt[iy * nx + ix],
361 ztm[iy * nx + ix] / npt[iy * nx + ix],
362 ttm[iy * nx + ix] / npt[iy * nx + ix],
363 h2otm[iy * nx + ix] / npt[iy * nx + ix],
364 lwcm[iy * nx + ix] / np[iy * nx + ix],
365 rwcm[iy * nx + ix] / np[iy * nx + ix],
366 iwcm[iy * nx + ix] / np[iy * nx + ix],
367 swcm[iy * nx + ix] / np[iy * nx + ix],
368 ccm[iy * nx + ix] / np[iy * nx + ix],
369 clm[iy * nx + ix] / np[iy * nx + ix],
370 pctm[iy * nx + ix] / np[iy * nx + ix],
371 pcbm[iy * nx + ix] / np[iy * nx + ix],
372 plclm[iy * nx + ix] / npc[iy * nx + ix],
373 plfcm[iy * nx + ix] / npc[iy * nx + ix],
374 pelm[iy * nx + ix] / npc[iy * nx + ix],
375 capem[iy * nx + ix] / npc[iy * nx + ix],
376 cinm[iy * nx + ix] / npc[iy * nx + ix],
377 rhm[iy * nx + ix] / np[iy * nx + ix],
378 rhicem[iy * nx + ix] / np[iy * nx + ix],
379 tdewm[iy * nx + ix] / np[iy * nx + ix],
380 ticem[iy * nx + ix] / np[iy * nx + ix],
381 tnatm[iy * nx + ix] / np[iy * nx + ix],
382 hno3m[iy * nx + ix] / np[iy * nx + ix],
383 ohm[iy * nx + ix] / np[iy * nx + ix],
384 h2o2m[iy * nx + ix] / np[iy * nx + ix],
385 ho2m[iy * nx + ix] / np[iy * nx + ix],
386 o1dm[iy * nx + ix] / np[iy * nx + ix],
387 pblm[iy * nx + ix] / np[iy * nx + ix],
388 o3cm[iy * nx + ix] / np[iy * nx + ix], np[iy * nx + ix],
389 npt[iy * nx + ix], npc[iy * nx + ix]);
390 }
391
392
393 fclose(out);
394
395
396 mptrac_free(NULL, NULL, clim, met, NULL, NULL, NULL, dd);
397 free(timem);
398 free(psm);
399 free(tsm);
400 free(zsm);
401 free(usm);
402 free(vsm);
403 free(essm);
404 free(nssm);
405 free(shfm);
406 free(lsmm);
407 free(sstm);
408 free(pblm);
409 free(ptm);
410 free(pm);
411 free(tm);
412 free(um);
413 free(vm);
414 free(wm);
415 free(h2om);
416 free(h2otm);
417 free(o3m);
418 free(hno3m);
419 free(ohm);
420 free(h2o2m);
421 free(ho2m);
422 free(o1dm);
423 free(tdewm);
424 free(ticem);
425 free(tnatm);
426 free(lwcm);
427 free(rwcm);
428 free(iwcm);
429 free(swcm);
430 free(ccm);
431 free(zm);
432 free(pvm);
433 free(ztm);
434 free(ttm);
435 free(pctm);
436 free(pcbm);
437 free(clm);
438 free(plclm);
439 free(plfcm);
440 free(pelm);
441 free(capem);
442 free(cinm);
443 free(o3cm);
444 free(rhm);
445 free(rhicem);
446 free(np);
447 free(npc);
448 free(npt);
449
450 return EXIT_SUCCESS;
451}
#define NX
Maximum number of longitudes.
#define NY
Maximum number of latitudes.
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.
void intpol_met_space_3d(const met_t *met, float array[EX][EY][EP], const double p, const double lon, const double lat, double *var, int *ci, double *cw, const int init)
Interpolates meteorological variables in 3D space.
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 INTPOL_INIT
Initialize arrays for interpolation.
#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 THETA(p, t)
Compute potential temperature.
#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 ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
#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.
float t[EX][EY][EP]
Temperature [K].
double lon[EX]
Longitudes [deg].
double lat[EY]
Latitudes [deg].
double p[EP]
Pressure levels [hPa].