53 static double *timem, ps, *psm, ts, *tsm, zs, *zsm, us, *usm, vs, *vsm, lsm,
54 *lsmm, sst, *sstm, pbl, *pblm, pt, *ptm, t, *pm, *tm, u, *um, v, *vm, w,
55 *wm, h2o, *h2om, h2ot, *h2otm, o3, *o3m, *hno3m, *ohm, *h2o2m, *ho2m,
56 *o1dm, *tdewm, *ticem, *tnatm, lwc, *lwcm, rwc, *rwcm, iwc, *iwcm, swc,
57 *swcm, cc, *ccm, z, *zm, pv, *pvm, zt, *ztm, tt, *ttm, pct, *pctm, pcb,
58 *pcbm, cl, *clm, plcl, *plclm, plfc, *plfcm, pel, *pelm, cape, *capem,
59 cin, *cinm, o3c, *o3cm, *rhm, *rhicem, ptop, pbot, t0, lon, lons[
NX], lat,
62 static int *np, *npc, *npt, nx, ny;
157 ALLOC(rhicem,
double,
168 ERRMSG(
"Give parameters: <ctl> <map.tab> <met0> [ <met1> ... ]");
171 read_ctl(argv[1], argc, argv, &ctl);
172 double p0 =
P(
scan_ctl(argv[1], argc, argv,
"MAP_Z0", -1,
"10", NULL));
173 double lon0 =
scan_ctl(argv[1], argc, argv,
"MAP_LON0", -1,
"-180", NULL);
174 double lon1 =
scan_ctl(argv[1], argc, argv,
"MAP_LON1", -1,
"180", NULL);
175 double dlon =
scan_ctl(argv[1], argc, argv,
"MAP_DLON", -1,
"-999", NULL);
176 double lat0 =
scan_ctl(argv[1], argc, argv,
"MAP_LAT0", -1,
"-90", NULL);
177 double lat1 =
scan_ctl(argv[1], argc, argv,
"MAP_LAT1", -1,
"90", NULL);
178 double dlat =
scan_ctl(argv[1], argc, argv,
"MAP_DLAT", -1,
"-999", NULL);
179 double theta =
scan_ctl(argv[1], argc, argv,
"MAP_THETA", -1,
"-999", NULL);
185 for (
int i = 3; i < argc; i++) {
188 if (!
read_met(argv[i], &ctl, clim, met))
193 dlon = fabs(met->
lon[1] - met->
lon[0]);
195 dlat = fabs(met->
lat[1] - met->
lat[0]);
196 if (lon0 < -360 && lon1 > 360) {
197 lon0 = gsl_stats_min(met->
lon, 1, (
size_t) met->
nx);
198 lon1 = gsl_stats_max(met->
lon, 1, (
size_t) met->
nx);
201 for (lon = lon0; lon <= lon1; lon += dlon) {
204 ERRMSG(
"Too many longitudes!");
206 if (lat0 < -90 && lat1 > 90) {
207 lat0 = gsl_stats_min(met->
lat, 1, (
size_t) met->
ny);
208 lat1 = gsl_stats_max(met->
lat, 1, (
size_t) met->
ny);
210 for (lat = lat0; lat <= lat1; lat += dlat) {
213 ERRMSG(
"Too many latitudes!");
217 for (
int ix = 0; ix < nx; ix++)
218 for (
int iy = 0; iy < ny; iy++) {
223 ptop = met->
p[met->
np - 1];
226 p0 = 0.5 * (ptop + pbot);
229 if (
THETA(p0, t0) > theta)
233 }
while (fabs(ptop - pbot) > 1e-5);
240 timem[iy * nx + ix] += met->
time;
241 zm[iy * nx + ix] += z;
242 pm[iy * nx + ix] += p0;
243 tm[iy * nx + ix] += t;
244 um[iy * nx + ix] += u;
245 vm[iy * nx + ix] += v;
246 wm[iy * nx + ix] += w;
247 pvm[iy * nx + ix] += pv;
248 h2om[iy * nx + ix] += h2o;
249 o3m[iy * nx + ix] += o3;
250 lwcm[iy * nx + ix] += lwc;
251 rwcm[iy * nx + ix] += rwc;
252 iwcm[iy * nx + ix] += iwc;
253 swcm[iy * nx + ix] += swc;
254 ccm[iy * nx + ix] += cc;
255 psm[iy * nx + ix] += ps;
256 tsm[iy * nx + ix] += ts;
257 zsm[iy * nx + ix] += zs;
258 usm[iy * nx + ix] += us;
259 vsm[iy * nx + ix] += vs;
260 lsmm[iy * nx + ix] += lsm;
261 sstm[iy * nx + ix] += sst;
262 pblm[iy * nx + ix] += pbl;
263 pctm[iy * nx + ix] += pct;
264 pcbm[iy * nx + ix] += pcb;
265 clm[iy * nx + ix] += cl;
266 if (isfinite(plfc) && isfinite(pel) && cape >= ctl.
conv_cape
268 plclm[iy * nx + ix] += plcl;
269 plfcm[iy * nx + ix] += plfc;
270 pelm[iy * nx + ix] += pel;
271 capem[iy * nx + ix] += cape;
272 cinm[iy * nx + ix] += cin;
276 ptm[iy * nx + ix] += pt;
277 ztm[iy * nx + ix] += zt;
278 ttm[iy * nx + ix] += tt;
279 h2otm[iy * nx + ix] += h2ot;
282 o3cm[iy * nx + ix] += o3c;
284 tnatm[iy * nx + ix] +=
288 clim_oh(&ctl, clim, met->
time, lons[ix], lats[iy], p0);
290 ho2m[iy * nx + ix] +=
clim_zm(&clim->
ho2, met->
time, lats[iy], p0);
291 o1dm[iy * nx + ix] +=
clim_zm(&clim->
o1d, met->
time, lats[iy], p0);
292 rhm[iy * nx + ix] +=
RH(p0, t, h2o);
293 rhicem[iy * nx + ix] +=
RHICE(p0, t, h2o);
294 tdewm[iy * nx + ix] +=
TDEW(p0, h2o);
295 ticem[iy * nx + ix] +=
TICE(p0, h2o);
301 LOG(1,
"Write meteorological data file: %s", argv[2]);
302 if (!(out = fopen(argv[2],
"w")))
303 ERRMSG(
"Cannot create file!");
309 for (
int iy = 0; iy < ny; iy++) {
311 for (
int ix = 0; ix < nx; ix++)
313 "%.2f %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
314 " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
315 " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %d %d %d\n",
316 timem[iy * nx + ix] / np[iy * nx + ix],
317 Z(pm[iy * nx + ix] / np[iy * nx + ix]), lons[ix], lats[iy],
318 pm[iy * nx + ix] / np[iy * nx + ix],
319 tm[iy * nx + ix] / np[iy * nx + ix],
320 um[iy * nx + ix] / np[iy * nx + ix],
321 vm[iy * nx + ix] / np[iy * nx + ix],
322 wm[iy * nx + ix] / np[iy * nx + ix],
323 h2om[iy * nx + ix] / np[iy * nx + ix],
324 o3m[iy * nx + ix] / np[iy * nx + ix],
325 zm[iy * nx + ix] / np[iy * nx + ix],
326 pvm[iy * nx + ix] / np[iy * nx + ix],
327 psm[iy * nx + ix] / np[iy * nx + ix],
328 tsm[iy * nx + ix] / np[iy * nx + ix],
329 zsm[iy * nx + ix] / np[iy * nx + ix],
330 usm[iy * nx + ix] / np[iy * nx + ix],
331 vsm[iy * nx + ix] / np[iy * nx + ix],
332 lsmm[iy * nx + ix] / np[iy * nx + ix],
333 sstm[iy * nx + ix] / np[iy * nx + ix],
334 ptm[iy * nx + ix] / npt[iy * nx + ix],
335 ztm[iy * nx + ix] / npt[iy * nx + ix],
336 ttm[iy * nx + ix] / npt[iy * nx + ix],
337 h2otm[iy * nx + ix] / npt[iy * nx + ix],
338 lwcm[iy * nx + ix] / np[iy * nx + ix],
339 rwcm[iy * nx + ix] / np[iy * nx + ix],
340 iwcm[iy * nx + ix] / np[iy * nx + ix],
341 swcm[iy * nx + ix] / np[iy * nx + ix],
342 ccm[iy * nx + ix] / np[iy * nx + ix],
343 clm[iy * nx + ix] / np[iy * nx + ix],
344 pctm[iy * nx + ix] / np[iy * nx + ix],
345 pcbm[iy * nx + ix] / np[iy * nx + ix],
346 plclm[iy * nx + ix] / npc[iy * nx + ix],
347 plfcm[iy * nx + ix] / npc[iy * nx + ix],
348 pelm[iy * nx + ix] / npc[iy * nx + ix],
349 capem[iy * nx + ix] / npc[iy * nx + ix],
350 cinm[iy * nx + ix] / npc[iy * nx + ix],
351 rhm[iy * nx + ix] / np[iy * nx + ix],
352 rhicem[iy * nx + ix] / np[iy * nx + ix],
353 tdewm[iy * nx + ix] / np[iy * nx + ix],
354 ticem[iy * nx + ix] / np[iy * nx + ix],
355 tnatm[iy * nx + ix] / np[iy * nx + ix],
356 hno3m[iy * nx + ix] / np[iy * nx + ix],
357 ohm[iy * nx + ix] / np[iy * nx + ix],
358 h2o2m[iy * nx + ix] / np[iy * nx + ix],
359 ho2m[iy * nx + ix] / np[iy * nx + ix],
360 o1dm[iy * nx + ix] / np[iy * nx + ix],
361 pblm[iy * nx + ix] / np[iy * nx + ix],
362 o3cm[iy * nx + ix] / np[iy * nx + ix], np[iy * nx + ix],
363 npt[iy * nx + ix], npc[iy * nx + ix]);
int main(int argc, char *argv[])
#define NX
Maximum number of longitudes.
#define NY
Maximum number of latitudes.
void read_clim(ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
double clim_zm(const clim_zm_t *zm, const double t, const double lat, const double p)
Interpolates monthly mean zonal mean climatological variables.
double nat_temperature(const double p, const double h2o, const double hno3)
Calculates the nitric acid trihydrate (NAT) temperature.
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 intpol_met_space_3d(met_t *met, float array[EX][EY][EP], double p, double lon, double lat, double *var, int *ci, double *cw, int init)
Interpolates meteorological variables in 3D space.
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.
int read_met(const char *filename, ctl_t *ctl, clim_t *clim, met_t *met)
Reads meteorological data from a file and populates the provided structures.
double scan_ctl(const char *filename, int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Scans a control file or command-line arguments for a specified variable.
MPTRAC library declarations.
#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 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].
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]
Longitude [deg].
double lat[EY]
Latitude [deg].
double p[EP]
Pressure levels [hPa].