MPTRAC
Macros | Functions
met_map.c File Reference

Extract map from meteorological data. More...

#include "mptrac.h"

Go to the source code of this file.

Macros

#define NX   EX
 Maximum number of longitudes. More...
 
#define NY   EY
 Maximum number of latitudes. More...
 

Functions

int main (int argc, char *argv[])
 

Detailed Description

Extract map from meteorological data.

Definition in file met_map.c.

Macro Definition Documentation

◆ NX

#define NX   EX

Maximum number of longitudes.

Definition at line 32 of file met_map.c.

◆ NY

#define NY   EY

Maximum number of latitudes.

Definition at line 35 of file met_map.c.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 41 of file met_map.c.

43 {
44
45 ctl_t ctl;
46
47 clim_t *clim;
48
49 met_t *met;
50
51 FILE *out;
52
53 static double *timem, ps, *psm, ts, *tsm, zs, *zsm, us, *usm, vs, *vsm,
54 ess, *essm, nss, *nssm, shf, *shfm, lsm, *lsmm, sst, *sstm, pbl, *pblm,
55 pt, *ptm, t, *pm, *tm, u, *um, v, *vm, w, *wm, h2o, *h2om, h2ot, *h2otm,
56 o3, *o3m, *hno3m, *ohm, *h2o2m, *ho2m, *o1dm, *tdewm, *ticem, *tnatm,
57 lwc, *lwcm, rwc, *rwcm, iwc, *iwcm, swc, *swcm, cc, *ccm, z, *zm,
58 pv, *pvm, zt, *ztm, tt, *ttm, pct, *pctm, pcb, *pcbm, cl, *clm,
59 plcl, *plclm, plfc, *plfcm, pel, *pelm, cape, *capem, cin, *cinm,
60 o3c, *o3cm, *rhm, *rhicem, ptop, pbot, t0, lon, lons[NX], lat, lats[NY];
61
62 static int *np, *npc, *npt, nx, ny;
63
64 /* Allocate... */
65 ALLOC(clim, clim_t, 1);
66 ALLOC(met, met_t, 1);
67 ALLOC(timem, double,
68 NX * NY);
69 ALLOC(psm, double,
70 NX * NY);
71 ALLOC(tsm, double,
72 NX * NY);
73 ALLOC(zsm, double,
74 NX * NY);
75 ALLOC(usm, double,
76 NX * NY);
77 ALLOC(vsm, double,
78 NX * NY);
79 ALLOC(essm, double,
80 NX * NY);
81 ALLOC(nssm, double,
82 NX * NY);
83 ALLOC(shfm, double,
84 NX * NY);
85 ALLOC(lsmm, double,
86 NX * NY);
87 ALLOC(sstm, double,
88 NX * NY);
89 ALLOC(pblm, double,
90 NX * NY);
91 ALLOC(ptm, double,
92 NX * NY);
93 ALLOC(pm, double,
94 NX * NY);
95 ALLOC(tm, double,
96 NX * NY);
97 ALLOC(um, double,
98 NX * NY);
99 ALLOC(vm, double,
100 NX * NY);
101 ALLOC(wm, double,
102 NX * NY);
103 ALLOC(h2om, double,
104 NX * NY);
105 ALLOC(h2otm, double,
106 NX * NY);
107 ALLOC(o3m, double,
108 NX * NY);
109 ALLOC(hno3m, double,
110 NX * NY);
111 ALLOC(ohm, double,
112 NX * NY);
113 ALLOC(h2o2m, double,
114 NX * NY);
115 ALLOC(ho2m, double,
116 NX * NY);
117 ALLOC(o1dm, double,
118 NX * NY);
119 ALLOC(tdewm, double,
120 NX * NY);
121 ALLOC(ticem, double,
122 NX * NY);
123 ALLOC(tnatm, double,
124 NX * NY);
125 ALLOC(lwcm, double,
126 NX * NY);
127 ALLOC(rwcm, double,
128 NX * NY);
129 ALLOC(iwcm, double,
130 NX * NY);
131 ALLOC(swcm, double,
132 NX * NY);
133 ALLOC(ccm, double,
134 NX * NY);
135 ALLOC(zm, double,
136 NX * NY);
137 ALLOC(pvm, double,
138 NX * NY);
139 ALLOC(ztm, double,
140 NX * NY);
141 ALLOC(ttm, double,
142 NX * NY);
143 ALLOC(pctm, double,
144 NX * NY);
145 ALLOC(pcbm, double,
146 NX * NY);
147 ALLOC(clm, double,
148 NX * NY);
149 ALLOC(plclm, double,
150 NX * NY);
151 ALLOC(plfcm, double,
152 NX * NY);
153 ALLOC(pelm, double,
154 NX * NY);
155 ALLOC(capem, double,
156 NX * NY);
157 ALLOC(cinm, double,
158 NX * NY);
159 ALLOC(o3cm, double,
160 NX * NY);
161 ALLOC(rhm, double,
162 NX * NY);
163 ALLOC(rhicem, double,
164 NX * NY);
165 ALLOC(np, int,
166 NX * NY);
167 ALLOC(npc, int,
168 NX * NY);
169 ALLOC(npt, int,
170 NX * NY);
171
172 /* Check arguments... */
173 if (argc < 4)
174 ERRMSG("Give parameters: <ctl> <map.tab> <met0> [ <met1> ... ]");
175
176 /* Read control parameters... */
177 read_ctl(argv[1], argc, argv, &ctl);
178 double p0 = P(scan_ctl(argv[1], argc, argv, "MAP_Z0", -1, "10", NULL));
179 double lon0 = scan_ctl(argv[1], argc, argv, "MAP_LON0", -1, "-180", NULL);
180 double lon1 = scan_ctl(argv[1], argc, argv, "MAP_LON1", -1, "180", NULL);
181 double dlon = scan_ctl(argv[1], argc, argv, "MAP_DLON", -1, "-999", NULL);
182 double lat0 = scan_ctl(argv[1], argc, argv, "MAP_LAT0", -1, "-90", NULL);
183 double lat1 = scan_ctl(argv[1], argc, argv, "MAP_LAT1", -1, "90", NULL);
184 double dlat = scan_ctl(argv[1], argc, argv, "MAP_DLAT", -1, "-999", NULL);
185 double theta = scan_ctl(argv[1], argc, argv, "MAP_THETA", -1, "-999", NULL);
186
187 /* Read climatological data... */
188 read_clim(&ctl, clim);
189
190 /* Loop over files... */
191 for (int i = 3; i < argc; i++) {
192
193 /* Read meteorological data... */
194 if (!read_met(argv[i], &ctl, clim, met))
195 continue;
196
197 /* Set horizontal grid... */
198 if (dlon <= 0)
199 dlon = fabs(met->lon[1] - met->lon[0]);
200 if (dlat <= 0)
201 dlat = fabs(met->lat[1] - met->lat[0]);
202 if (lon0 < -360 && lon1 > 360) {
203 lon0 = gsl_stats_min(met->lon, 1, (size_t) met->nx);
204 lon1 = gsl_stats_max(met->lon, 1, (size_t) met->nx);
205 }
206 nx = ny = 0;
207 for (lon = lon0; lon <= lon1; lon += dlon) {
208 lons[nx] = lon;
209 if ((++nx) >= NX)
210 ERRMSG("Too many longitudes!");
211 }
212 if (lat0 < -90 && lat1 > 90) {
213 lat0 = gsl_stats_min(met->lat, 1, (size_t) met->ny);
214 lat1 = gsl_stats_max(met->lat, 1, (size_t) met->ny);
215 }
216 for (lat = lat0; lat <= lat1; lat += dlat) {
217 lats[ny] = lat;
218 if ((++ny) >= NY)
219 ERRMSG("Too many latitudes!");
220 }
221
222 /* Average... */
223 for (int ix = 0; ix < nx; ix++)
224 for (int iy = 0; iy < ny; iy++) {
225
226 /* Find pressure level for given theta level... */
227 INTPOL_INIT;
228 if (theta > 0) {
229 ptop = met->p[met->np - 1];
230 pbot = met->p[0];
231 do {
232 p0 = 0.5 * (ptop + pbot);
233 intpol_met_space_3d(met, met->t, p0, lons[ix], lats[iy],
234 &t0, ci, cw, 1);
235 if (THETA(p0, t0) > theta)
236 ptop = p0;
237 else
238 pbot = p0;
239 } while (fabs(ptop - pbot) > 1e-5);
240 }
241
242 /* Interpolate meteo data... */
243 INTPOL_SPACE_ALL(p0, lons[ix], lats[iy]);
244
245 /* Averaging... */
246 timem[iy * nx + ix] += met->time;
247 zm[iy * nx + ix] += z;
248 pm[iy * nx + ix] += p0;
249 tm[iy * nx + ix] += t;
250 um[iy * nx + ix] += u;
251 vm[iy * nx + ix] += v;
252 wm[iy * nx + ix] += w;
253 pvm[iy * nx + ix] += pv;
254 h2om[iy * nx + ix] += h2o;
255 o3m[iy * nx + ix] += o3;
256 lwcm[iy * nx + ix] += lwc;
257 rwcm[iy * nx + ix] += rwc;
258 iwcm[iy * nx + ix] += iwc;
259 swcm[iy * nx + ix] += swc;
260 ccm[iy * nx + ix] += cc;
261 psm[iy * nx + ix] += ps;
262 tsm[iy * nx + ix] += ts;
263 zsm[iy * nx + ix] += zs;
264 usm[iy * nx + ix] += us;
265 vsm[iy * nx + ix] += vs;
266 essm[iy * nx + ix] += ess;
267 nssm[iy * nx + ix] += nss;
268 shfm[iy * nx + ix] += shf;
269 lsmm[iy * nx + ix] += lsm;
270 sstm[iy * nx + ix] += sst;
271 pblm[iy * nx + ix] += pbl;
272 pctm[iy * nx + ix] += pct;
273 pcbm[iy * nx + ix] += pcb;
274 clm[iy * nx + ix] += cl;
275 if (isfinite(plfc) && isfinite(pel) && cape >= ctl.conv_cape
276 && (ctl.conv_cin <= 0 || cin < ctl.conv_cin)) {
277 plclm[iy * nx + ix] += plcl;
278 plfcm[iy * nx + ix] += plfc;
279 pelm[iy * nx + ix] += pel;
280 capem[iy * nx + ix] += cape;
281 cinm[iy * nx + ix] += cin;
282 npc[iy * nx + ix]++;
283 }
284 if (isfinite(pt)) {
285 ptm[iy * nx + ix] += pt;
286 ztm[iy * nx + ix] += zt;
287 ttm[iy * nx + ix] += tt;
288 h2otm[iy * nx + ix] += h2ot;
289 npt[iy * nx + ix]++;
290 }
291 o3cm[iy * nx + ix] += o3c;
292 hno3m[iy * nx + ix] += clim_zm(&clim->hno3, met->time, lats[iy], p0);
293 tnatm[iy * nx + ix] +=
294 nat_temperature(p0, h2o,
295 clim_zm(&clim->hno3, met->time, lats[iy], p0));
296 ohm[iy * nx + ix] +=
297 clim_oh(&ctl, clim, met->time, lons[ix], lats[iy], p0);
298 h2o2m[iy * nx + ix] += clim_zm(&clim->h2o2, met->time, lats[iy], p0);
299 ho2m[iy * nx + ix] += clim_zm(&clim->ho2, met->time, lats[iy], p0);
300 o1dm[iy * nx + ix] += clim_zm(&clim->o1d, met->time, lats[iy], p0);
301 rhm[iy * nx + ix] += RH(p0, t, h2o);
302 rhicem[iy * nx + ix] += RHICE(p0, t, h2o);
303 tdewm[iy * nx + ix] += TDEW(p0, h2o);
304 ticem[iy * nx + ix] += TICE(p0, h2o);
305 np[iy * nx + ix]++;
306 }
307 }
308
309 /* Create output file... */
310 LOG(1, "Write meteorological data file: %s", argv[2]);
311 if (!(out = fopen(argv[2], "w")))
312 ERRMSG("Cannot create file!");
313
314 /* Write header... */
315 MET_HEADER;
316
317 /* Write data... */
318 for (int iy = 0; iy < ny; iy++) {
319 fprintf(out, "\n");
320 for (int ix = 0; ix < nx; ix++)
321 fprintf(out,
322 "%.2f %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
323 " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
324 " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %d %d %d\n",
325 timem[iy * nx + ix] / np[iy * nx + ix],
326 Z(pm[iy * nx + ix] / np[iy * nx + ix]), lons[ix], lats[iy],
327 pm[iy * nx + ix] / np[iy * nx + ix],
328 tm[iy * nx + ix] / np[iy * nx + ix],
329 um[iy * nx + ix] / np[iy * nx + ix],
330 vm[iy * nx + ix] / np[iy * nx + ix],
331 wm[iy * nx + ix] / np[iy * nx + ix],
332 h2om[iy * nx + ix] / np[iy * nx + ix],
333 o3m[iy * nx + ix] / np[iy * nx + ix],
334 zm[iy * nx + ix] / np[iy * nx + ix],
335 pvm[iy * nx + ix] / np[iy * nx + ix],
336 psm[iy * nx + ix] / np[iy * nx + ix],
337 tsm[iy * nx + ix] / np[iy * nx + ix],
338 zsm[iy * nx + ix] / np[iy * nx + ix],
339 usm[iy * nx + ix] / np[iy * nx + ix],
340 vsm[iy * nx + ix] / np[iy * nx + ix],
341 essm[iy * nx + ix] / np[iy * nx + ix],
342 nssm[iy * nx + ix] / np[iy * nx + ix],
343 shfm[iy * nx + ix] / np[iy * nx + ix],
344 lsmm[iy * nx + ix] / np[iy * nx + ix],
345 sstm[iy * nx + ix] / np[iy * nx + ix],
346 ptm[iy * nx + ix] / npt[iy * nx + ix],
347 ztm[iy * nx + ix] / npt[iy * nx + ix],
348 ttm[iy * nx + ix] / npt[iy * nx + ix],
349 h2otm[iy * nx + ix] / npt[iy * nx + ix],
350 lwcm[iy * nx + ix] / np[iy * nx + ix],
351 rwcm[iy * nx + ix] / np[iy * nx + ix],
352 iwcm[iy * nx + ix] / np[iy * nx + ix],
353 swcm[iy * nx + ix] / np[iy * nx + ix],
354 ccm[iy * nx + ix] / np[iy * nx + ix],
355 clm[iy * nx + ix] / np[iy * nx + ix],
356 pctm[iy * nx + ix] / np[iy * nx + ix],
357 pcbm[iy * nx + ix] / np[iy * nx + ix],
358 plclm[iy * nx + ix] / npc[iy * nx + ix],
359 plfcm[iy * nx + ix] / npc[iy * nx + ix],
360 pelm[iy * nx + ix] / npc[iy * nx + ix],
361 capem[iy * nx + ix] / npc[iy * nx + ix],
362 cinm[iy * nx + ix] / npc[iy * nx + ix],
363 rhm[iy * nx + ix] / np[iy * nx + ix],
364 rhicem[iy * nx + ix] / np[iy * nx + ix],
365 tdewm[iy * nx + ix] / np[iy * nx + ix],
366 ticem[iy * nx + ix] / np[iy * nx + ix],
367 tnatm[iy * nx + ix] / np[iy * nx + ix],
368 hno3m[iy * nx + ix] / np[iy * nx + ix],
369 ohm[iy * nx + ix] / np[iy * nx + ix],
370 h2o2m[iy * nx + ix] / np[iy * nx + ix],
371 ho2m[iy * nx + ix] / np[iy * nx + ix],
372 o1dm[iy * nx + ix] / np[iy * nx + ix],
373 pblm[iy * nx + ix] / np[iy * nx + ix],
374 o3cm[iy * nx + ix] / np[iy * nx + ix], np[iy * nx + ix],
375 npt[iy * nx + ix], npc[iy * nx + ix]);
376 }
377
378 /* Close file... */
379 fclose(out);
380
381 /* Free... */
382 free(clim);
383 free(met);
384 free(timem);
385 free(psm);
386 free(tsm);
387 free(zsm);
388 free(usm);
389 free(vsm);
390 free(essm);
391 free(nssm);
392 free(shfm);
393 free(lsmm);
394 free(sstm);
395 free(pblm);
396 free(ptm);
397 free(pm);
398 free(tm);
399 free(um);
400 free(vm);
401 free(wm);
402 free(h2om);
403 free(h2otm);
404 free(o3m);
405 free(hno3m);
406 free(ohm);
407 free(h2o2m);
408 free(ho2m);
409 free(o1dm);
410 free(tdewm);
411 free(ticem);
412 free(tnatm);
413 free(lwcm);
414 free(rwcm);
415 free(iwcm);
416 free(swcm);
417 free(ccm);
418 free(zm);
419 free(pvm);
420 free(ztm);
421 free(ttm);
422 free(pctm);
423 free(pcbm);
424 free(clm);
425 free(plclm);
426 free(plfcm);
427 free(pelm);
428 free(capem);
429 free(cinm);
430 free(o3cm);
431 free(rhm);
432 free(rhicem);
433 free(np);
434 free(npc);
435 free(npt);
436
437 return EXIT_SUCCESS;
438}
NX
#define NX
Maximum number of longitudes.
Definition: met_map.c:32
NY
#define NY
Maximum number of latitudes.
Definition: met_map.c:35
read_met
int read_met(const char *filename, const ctl_t *ctl, const clim_t *clim, met_t *met)
Reads meteorological data from a file, supporting multiple formats and MPI broadcasting.
Definition: mptrac.c:6017
clim_zm
double clim_zm(const clim_zm_t *zm, const double t, const double lat, const double p)
Interpolates monthly mean zonal mean climatological variables.
Definition: mptrac.c:401
nat_temperature
double nat_temperature(const double p, const double h2o, const double hno3)
Calculates the nitric acid trihydrate (NAT) temperature.
Definition: mptrac.c:4520
intpol_met_space_3d
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.
Definition: mptrac.c:1459
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:8503
clim_oh
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...
Definition: mptrac.c:89
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:5156
read_clim
void read_clim(const ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
Definition: mptrac.c:4824
INTPOL_SPACE_ALL
#define INTPOL_SPACE_ALL(p, lon, lat)
Interpolate multiple meteorological variables in space.
Definition: mptrac.h:732
INTPOL_INIT
#define INTPOL_INIT
Initialize arrays for interpolation.
Definition: mptrac.h:682
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:1916
Z
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:1741
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1304
THETA
#define THETA(p, t)
Compute potential temperature.
Definition: mptrac.h:1644
MET_HEADER
#define MET_HEADER
Write header for meteorological data file.
Definition: mptrac.h:900
TICE
#define TICE(p, h2o)
Calculate frost point temperature (WMO, 2018).
Definition: mptrac.h:1620
RHICE
#define RHICE(p, t, h2o)
Compute relative humidity over ice.
Definition: mptrac.h:1456
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:349
RH
#define RH(p, t, h2o)
Compute relative humidity over water.
Definition: mptrac.h:1426
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:1846
TDEW
#define TDEW(p, h2o)
Calculate dew point temperature.
Definition: mptrac.h:1595
clim_t
Climatological data.
Definition: mptrac.h:3315
clim_t::ho2
clim_zm_t ho2
HO2 zonal means.
Definition: mptrac.h:3345
clim_t::hno3
clim_zm_t hno3
HNO3 zonal means.
Definition: mptrac.h:3336
clim_t::o1d
clim_zm_t o1d
O(1D) zonal means.
Definition: mptrac.h:3348
clim_t::h2o2
clim_zm_t h2o2
H2O2 zonal means.
Definition: mptrac.h:3342
ctl_t
Control parameters.
Definition: mptrac.h:2170
ctl_t::conv_cape
double conv_cape
CAPE threshold for convection module [J/kg].
Definition: mptrac.h:2693
ctl_t::conv_cin
double conv_cin
CIN threshold for convection module [J/kg].
Definition: mptrac.h:2696
met_t
Meteo data structure.
Definition: mptrac.h:3374
met_t::nx
int nx
Number of longitudes.
Definition: mptrac.h:3380
met_t::ny
int ny
Number of latitudes.
Definition: mptrac.h:3383
met_t::np
int np
Number of pressure levels.
Definition: mptrac.h:3386
met_t::t
float t[EX][EY][EP]
Temperature [K].
Definition: mptrac.h:3479
met_t::lon
double lon[EX]
Longitude [deg].
Definition: mptrac.h:3392
met_t::time
double time
Time [s].
Definition: mptrac.h:3377
met_t::lat
double lat[EY]
Latitude [deg].
Definition: mptrac.h:3395
met_t::p
double p[EP]
Pressure levels [hPa].
Definition: mptrac.h:3398
Here is the call graph for this function:
Generated on Wed Jan 22 2025 18:02:19 for MPTRAC by doxygen 1.9.4