MPTRAC
atm_init.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-2026 Forschungszentrum Juelich GmbH
18*/
19
25#include "mptrac.h"
26
27/* ------------------------------------------------------------
28 Functions...
29 ------------------------------------------------------------ */
30
32void usage(
33 void);
34
35/* ------------------------------------------------------------
36 Main...
37 ------------------------------------------------------------ */
38
39int main(
40 int argc,
41 char *argv[]) {
42
43 ctl_t ctl;
44
45 clim_t *clim;
46
47 atm_t *atm;
48
49 met_t *met0, *met1;
50
51 dd_t *dd;
52
53 /* Print usage information... */
54 USAGE;
55
56 /* Check arguments... */
57 if (argc < 3)
58 ERRMSG("Missing or invalid command-line arguments.\n\n"
59 "Usage: atm_init <ctl> <atm_out> [KEY VALUE ...]\n\n"
60 "Use -h for full help.");
61
62 /* Allocate... */
63 ALLOC(clim, clim_t, 1);
64 ALLOC(atm, atm_t, 1);
65 ALLOC(met0, met_t, 1);
66 ALLOC(met1, met_t, 1);
67 ALLOC(dd, dd_t, 1);
68
69 /* Read control parameters... */
70 mptrac_read_ctl(argv[1], argc, argv, &ctl);
71 const int ens =
72 (int) scan_ctl(argv[1], argc, argv, "INIT_ENS", -1, "0", NULL);
73 const double t0 = scan_ctl(argv[1], argc, argv, "INIT_T0", -1, "0", NULL);
74 const double t1 = scan_ctl(argv[1], argc, argv, "INIT_T1", -1, "0", NULL);
75 const double dt = scan_ctl(argv[1], argc, argv, "INIT_DT", -1, "1", NULL);
76 const double z0 = scan_ctl(argv[1], argc, argv, "INIT_Z0", -1, "0", NULL);
77 const double z1 = scan_ctl(argv[1], argc, argv, "INIT_Z1", -1, "0", NULL);
78 const double dz = scan_ctl(argv[1], argc, argv, "INIT_DZ", -1, "1", NULL);
79 const double relz0 =
80 scan_ctl(argv[1], argc, argv, "INIT_RELZ0", -1, "-999", NULL);
81 const double relz1 =
82 scan_ctl(argv[1], argc, argv, "INIT_RELZ1", -1, "-999", NULL);
83 const double drelz =
84 scan_ctl(argv[1], argc, argv, "INIT_DRELZ", -1, "1", NULL);
85 const int random_relz =
86 (int) scan_ctl(argv[1], argc, argv, "INIT_RANDOM_RELZ", -1, "0", NULL);
87 const int use_relz = (relz0 > -900.0 || relz1 > -900.0);
88 const double lon0 =
89 scan_ctl(argv[1], argc, argv, "INIT_LON0", -1, "0", NULL);
90 const double lon1 =
91 scan_ctl(argv[1], argc, argv, "INIT_LON1", -1, "0", NULL);
92 const double dlon =
93 scan_ctl(argv[1], argc, argv, "INIT_DLON", -1, "1", NULL);
94 const double lat0 =
95 scan_ctl(argv[1], argc, argv, "INIT_LAT0", -1, "0", NULL);
96 const double lat1 =
97 scan_ctl(argv[1], argc, argv, "INIT_LAT1", -1, "0", NULL);
98 const double dlat =
99 scan_ctl(argv[1], argc, argv, "INIT_DLAT", -1, "1", NULL);
100 const double st = scan_ctl(argv[1], argc, argv, "INIT_ST", -1, "0", NULL);
101 const double sz = scan_ctl(argv[1], argc, argv, "INIT_SZ", -1, "0", NULL);
102 const double slon =
103 scan_ctl(argv[1], argc, argv, "INIT_SLON", -1, "0", NULL);
104 const double slat =
105 scan_ctl(argv[1], argc, argv, "INIT_SLAT", -1, "0", NULL);
106 const double sx = scan_ctl(argv[1], argc, argv, "INIT_SX", -1, "0", NULL);
107 const double ut = scan_ctl(argv[1], argc, argv, "INIT_UT", -1, "0", NULL);
108 const double uz = scan_ctl(argv[1], argc, argv, "INIT_UZ", -1, "0", NULL);
109 const double ulon =
110 scan_ctl(argv[1], argc, argv, "INIT_ULON", -1, "0", NULL);
111 const double ulat =
112 scan_ctl(argv[1], argc, argv, "INIT_ULAT", -1, "0", NULL);
113 const int even =
114 (int) scan_ctl(argv[1], argc, argv, "INIT_EVENLY", -1, "0", NULL);
115 const int col_mass =
116 (int) scan_ctl(argv[1], argc, argv, "INIT_COL_MASS", -1, "0", NULL);
117 const int well_mixed =
118 (int) scan_ctl(argv[1], argc, argv, "INIT_WELL_MIXED", -1, "0", NULL);
119 const int rep =
120 (int) scan_ctl(argv[1], argc, argv, "INIT_REP", -1, "1", NULL);
121 const int init_np =
122 (int) scan_ctl(argv[1], argc, argv, "INIT_NP", -1, "0", NULL);
123 const double m = scan_ctl(argv[1], argc, argv, "INIT_MASS", -1, "0", NULL);
124 const double vmr = scan_ctl(argv[1], argc, argv, "INIT_VMR", -1, "0", NULL);
125 const double bellrad =
126 scan_ctl(argv[1], argc, argv, "INIT_BELLRAD", -1, "0", NULL);
127 const int idx_offset =
128 (int) scan_ctl(argv[1], argc, argv, "INIT_IDX_OFFSET", -1, "0", NULL);
129
130 /* Check arguments... */
131 if (ctl.met_coord_type != 0 && even)
132 ERRMSG("INIT_EVENLY is only supported for lat/lon grids");
133 if (ctl.met_coord_type != 0 && bellrad > 0)
134 ERRMSG("INIT_BELLRAD is only supported for lat/lon grids");
135 if (rep <= 0)
136 ERRMSG("INIT_REP must be positive!");
137 if (init_np < 0 || init_np > NP)
138 ERRMSG("INIT_NP must be between 0 and NP!");
139 if (well_mixed && ctl.met_coord_type != 0)
140 ERRMSG("INIT_WELL_MIXED is only supported for lat/lon grids!");
141 if (well_mixed && !even)
142 ERRMSG("INIT_WELL_MIXED requires INIT_EVENLY=1!");
143 if (well_mixed && !col_mass)
144 ERRMSG("INIT_WELL_MIXED requires INIT_COL_MASS=1!");
145 if (init_np > 0 && !well_mixed)
146 ERRMSG("INIT_NP requires INIT_WELL_MIXED!");
147 if (col_mass && (t1 != t0 || st > 0 || ut > 0))
148 ERRMSG
149 ("INIT_COL_MASS/INIT_WELL_MIXED require INIT_T0 == INIT_T1 and no time scattering!");
150 if (use_relz && (z0 != 0.0 || z1 != 0.0))
151 ERRMSG("INIT_RELZ0/INIT_RELZ1 cannot be combined with INIT_Z0/INIT_Z1!");
152 if (use_relz && (relz0 <= -900.0 || relz1 <= -900.0))
153 ERRMSG("INIT_RELZ0 and INIT_RELZ1 must both be set!");
154 if (use_relz && (t1 != t0 || st > 0 || ut > 0))
155 ERRMSG
156 ("INIT_RELZ0/INIT_RELZ1 require INIT_T0 == INIT_T1 and no time scattering!");
157 if (use_relz && random_relz && relz1 < relz0)
158 ERRMSG("INIT_RELZ1 must be greater than or equal to INIT_RELZ0!");
159 if (use_relz && !random_relz && drelz <= 0)
160 ERRMSG("INIT_DRELZ must be positive!");
161 if (!use_relz && dz <= 0)
162 ERRMSG("INIT_DZ must be positive!");
163
164 /* Read climatological data... */
165 mptrac_read_clim(&ctl, clim);
166
167 /* Initialize random number generator... */
168 gsl_rng_env_setup();
169 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
170
171 /* Get meteorological data for column-mass, well-mixed, and relative-PBL initialization... */
172 double ptop = 0, dpcolmax = 0;
173 if (col_mass || use_relz) {
174 mptrac_get_met(&ctl, clim, t0, &met0, &met1, dd);
175 ptop = gsl_stats_min(met0->p, 1, (size_t) met0->np);
176 double psmax = -1;
177 for (int ix = 0; ix < met0->nx; ix++)
178 for (int iy = 0; iy < met0->ny; iy++) {
179 psmax = MAX(psmax, met0->ps[ix][iy]);
180 psmax = MAX(psmax, met1->ps[ix][iy]);
181 }
182 dpcolmax = psmax - ptop;
183 }
184
185 /* Create grid... */
186 int ntry = 0;
187 const int ntrymax = init_np > 0 ? 1000000 * init_np : 0;
188 do {
189 const double vert0 = use_relz ? relz0 : z0;
190 const double vert1 = use_relz ? (random_relz ? relz0 : relz1) : z1;
191 const double dvert = use_relz ? (random_relz ? 1.0 : drelz) : dz;
192 for (double t = t0; t <= t1 && !(init_np > 0 && atm->np >= init_np);
193 t += dt)
194 for (double vert = vert0;
195 vert <= vert1 && !(init_np > 0 && atm->np >= init_np);
196 vert += dvert)
197 for (double lon = lon0;
198 lon <= lon1 && !(init_np > 0 && atm->np >= init_np); lon += dlon)
199 for (double lat = lat0;
200 lat <= lat1 && !(init_np > 0 && atm->np >= init_np);
201 lat += dlat)
202 for (int irep = 0;
203 irep < rep && !(init_np > 0 && atm->np >= init_np); irep++) {
204
205 /* Check number of attempts... */
206 if (init_np > 0 && ++ntry > ntrymax)
207 ERRMSG
208 ("INIT_NP could not be reached by INIT_WELL_MIXED sampling!");
209
210 /* Set time... */
211 double rg = gsl_ran_gaussian_ziggurat(rng, st / 2.3548);
212 double ru = ut * (gsl_rng_uniform(rng) - 0.5);
213 atm->time[atm->np] = (t + rg + ru);
214
215 if (!use_relz) {
216 /* Keep the legacy RNG order for absolute-height initialization. */
217 rg = gsl_ran_gaussian_ziggurat(rng, sz / 2.3548);
218 ru = uz * (gsl_rng_uniform(rng) - 0.5);
219 atm->p[atm->np] = P(vert + rg + ru);
220 }
221
222 /* Set horizontal position... */
223 rg = gsl_ran_gaussian_ziggurat(rng, slon / 2.3548);
224 const double sx_coord =
225 ctl.met_coord_type == 0 ? DX2DEG(sx, lat) : sx;
226 double rx = gsl_ran_gaussian_ziggurat(rng, sx_coord / 2.3548);
227 ru = ulon * (gsl_rng_uniform(rng) - 0.5);
228 atm->lon[atm->np] = (lon + rg + rx + ru);
229
230 /* Set ensemble index... */
231 if (ctl.qnt_ens >= 0)
232 atm->q[ctl.qnt_ens][atm->np] = ens;
233
234 /* Apply cosine-latitude weighting... */
235 do {
236 rg = gsl_ran_gaussian_ziggurat(rng, slat / 2.3548);
237 const double sy_coord =
238 ctl.met_coord_type == 0 ? DY2DEG(sx) : sx;
239 rx = gsl_ran_gaussian_ziggurat(rng, sy_coord / 2.3548);
240 ru = ulat * (gsl_rng_uniform(rng) - 0.5);
241 atm->lat[atm->np] = (lat + rg + rx + ru);
242 } while (even && gsl_rng_uniform(rng) >
243 fabs(cos(DEG2RAD(atm->lat[atm->np]))));
244
245 if (use_relz) {
246 /* Relative-PBL initialization needs the final horizontal position first. */
247 rg = gsl_ran_gaussian_ziggurat(rng, sz / 2.3548);
248 ru = uz * (gsl_rng_uniform(rng) - 0.5);
249 const double relz_sample =
250 random_relz ? relz0 + gsl_rng_uniform(rng) * (relz1 - relz0)
251 : vert;
252 INTPOL_INIT;
253 double ps, pbl;
254 intpol_met_time_2d(met0, met0->ps, met1, met1->ps,
255 atm->time[atm->np], atm->lon[atm->np],
256 atm->lat[atm->np], &ps, ci, cw, 1);
257 intpol_met_time_2d(met0, met0->pbl, met1, met1->pbl,
258 atm->time[atm->np], atm->lon[atm->np],
259 atm->lat[atm->np], &pbl, ci, cw, 1);
260 if (ctl.qnt_ps >= 0)
261 atm->q[ctl.qnt_ps][atm->np] = ps;
262 if (ctl.qnt_pbl >= 0)
263 atm->q[ctl.qnt_pbl][atm->np] = pbl;
264 const double zsurf = Z(ps);
265 const double zpbl = Z(pbl);
266 const double depth = zpbl - zsurf;
267 if (depth <= 0)
268 continue;
269 atm->p[atm->np] = P(zsurf + (relz_sample + rg + ru) * depth);
270 }
271
272 /* Apply cosine bell (Williamson et al., 1992)... */
273 if (bellrad > 0) {
274 double x0[3], x1[3];
275 geo2cart(0.0, 0.5 * (lon0 + lon1), 0.5 * (lat0 + lat1), x0);
276 geo2cart(0.0, atm->lon[atm->np], atm->lat[atm->np], x1);
277 const double rad =
278 RE * acos(DOTP(x0, x1) / sqrt(DOTP(x0, x0)) /
279 sqrt(DOTP(x1, x1)));
280 if (rad > bellrad)
281 continue;
282 if (ctl.qnt_m >= 0)
283 atm->q[ctl.qnt_m][atm->np] =
284 0.5 * (1. + cos(M_PI * rad / bellrad));
285 if (ctl.qnt_vmr >= 0)
286 atm->q[ctl.qnt_vmr][atm->np] =
287 0.5 * (1. + cos(M_PI * rad / bellrad));
288 }
289
290 /* Apply column-mass and globally well-mixed air-mass weighting... */
291 if (col_mass) {
292 INTPOL_INIT;
293 double ps;
294 intpol_met_time_2d(met0, met0->ps, met1, met1->ps,
295 atm->time[atm->np], atm->lon[atm->np],
296 atm->lat[atm->np], &ps, ci, cw, 1);
297 const double dpcol = ps - ptop;
298 if (gsl_rng_uniform(rng) > dpcol / dpcolmax)
299 continue;
300 if (well_mixed)
301 atm->p[atm->np] = ptop + dpcol * gsl_rng_uniform(rng);
302 }
303
304 /* Set particle counter... */
305 if ((++atm->np) > NP)
306 ERRMSG("Too many particles!");
307 }
308
309 } while (init_np > 0 && atm->np < init_np);
310
311 /* Check number of air parcels... */
312 if (atm->np <= 0)
313 ERRMSG("Did not create any air parcels!");
314
315 /* Initialize mass... */
316 if (ctl.qnt_m >= 0 && bellrad <= 0)
317 for (int ip = 0; ip < atm->np; ip++)
318 atm->q[ctl.qnt_m][ip] = m / atm->np;
319
320 /* Initialize volume mixing ratio... */
321 if (ctl.qnt_vmr >= 0 && bellrad <= 0)
322 for (int ip = 0; ip < atm->np; ip++)
323 atm->q[ctl.qnt_vmr][ip] = vmr;
324
325 /* Initialize air parcel index... */
326 if (ctl.qnt_idx >= 0)
327 for (int ip = 0; ip < atm->np; ip++)
328 atm->q[ctl.qnt_idx][ip] = idx_offset + ip;
329
330 /* Initialize age of air... */
331 if (ctl.qnt_aoa >= 0)
332 for (int ip = 0; ip < atm->np; ip++)
333 atm->q[ctl.qnt_aoa][ip] = atm->time[ip];
334
335 /* Save data... */
336 mptrac_write_atm(argv[2], &ctl, atm, 0);
337
338 /* Free... */
339 gsl_rng_free(rng);
340 free(clim);
341 free(atm);
342 free(met0);
343 free(met1);
344 free(dd);
345
346 return EXIT_SUCCESS;
347}
348
349/*****************************************************************************/
350
352void usage(
353 void) {
354
355 printf("\nMPTRAC atm_init tool.\n\n");
356 printf
357 ("Create an atmospheric data file with initial air parcel positions.\n");
358 printf("\n");
359 printf("Usage:\n");
360 printf(" atm_init <ctl> <atm_out> [KEY VALUE ...]\n");
361 printf("\n");
362 printf("Arguments:\n");
363 printf(" <ctl> Control file.\n");
364 printf(" <atm_out> Atmospheric output file.\n");
365 printf(" [KEY VALUE] Optional control parameters.\n");
366 printf("\nFurther information:\n");
367 printf(" Manual: https://slcs-jsc.github.io/mptrac/\n");
368}
main
int main(int argc, char *argv[])
Definition: atm_init.c:39
usage
void usage(void)
Print command-line help.
Definition: atm_init.c:352
mptrac_write_atm
void mptrac_write_atm(const char *filename, const ctl_t *ctl, const atm_t *atm, const double t)
Writes air parcel data to a file in various formats.
Definition: mptrac.c:7806
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:12110
mptrac_get_met
void mptrac_get_met(ctl_t *ctl, clim_t *clim, const double t, met_t **met0, met_t **met1, dd_t *dd)
Retrieves meteorological data for the specified time.
Definition: mptrac.c:6187
mptrac_read_clim
void mptrac_read_clim(const ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
Definition: mptrac.c:6406
intpol_met_time_2d
void intpol_met_time_2d(const met_t *met0, float array0[EX][EY], const met_t *met1, float array1[EX][EY], const double ts, const double lon, const double lat, double *var, int *ci, double *cw, const int init)
Interpolates meteorological data in 2D space and time.
Definition: mptrac.c:3141
mptrac_read_ctl
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.
Definition: mptrac.c:6466
geo2cart
void geo2cart(const double z, const double lon, const double lat, double *x)
Converts geographic coordinates (longitude, latitude, altitude) to Cartesian coordinates.
Definition: mptrac.c:2601
mptrac.h
MPTRAC library declarations.
RE
#define RE
Mean radius of Earth [km].
Definition: mptrac.h:315
DOTP
#define DOTP(a, b)
Calculate the dot product of two vectors.
Definition: mptrac.h:826
INTPOL_INIT
#define INTPOL_INIT
Initialize arrays for interpolation.
Definition: mptrac.h:964
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:2196
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: mptrac.h:2003
Z
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:2033
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1574
DX2DEG
#define DX2DEG(dx, lat)
Convert a distance in kilometers to degrees longitude at a given latitude.
Definition: mptrac.h:694
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:457
DEG2RAD
#define DEG2RAD(deg)
Converts degrees to radians.
Definition: mptrac.h:647
NP
#define NP
Maximum number of atmospheric data points.
Definition: mptrac.h:359
DY2DEG
#define DY2DEG(dy)
Convert a distance in kilometers to degrees latitude.
Definition: mptrac.h:712
MAX
#define MAX(a, b)
Macro to determine the maximum of two values.
Definition: mptrac.h:1168
atm_t
Air parcel data.
Definition: mptrac.h:3341
atm_t::time
double time[NP]
Time [s].
Definition: mptrac.h:3347
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3356
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3353
atm_t::np
int np
Number of air parcels.
Definition: mptrac.h:3344
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3359
atm_t::p
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3350
clim_t
Climatological data.
Definition: mptrac.h:3536
ctl_t
Control parameters.
Definition: mptrac.h:2284
ctl_t::qnt_m
int qnt_m
Quantity array index for mass.
Definition: mptrac.h:2315
ctl_t::qnt_aoa
int qnt_aoa
Quantity array index for age of air.
Definition: mptrac.h:2585
ctl_t::qnt_vmr
int qnt_vmr
Quantity array index for volume mixing ratio.
Definition: mptrac.h:2318
ctl_t::qnt_pbl
int qnt_pbl
Quantity array index for boundary layer pressure.
Definition: mptrac.h:2357
ctl_t::qnt_ps
int qnt_ps
Quantity array index for surface pressure.
Definition: mptrac.h:2327
ctl_t::qnt_ens
int qnt_ens
Quantity array index for ensemble IDs.
Definition: mptrac.h:2309
ctl_t::met_coord_type
int met_coord_type
Type of coordinates for meteo data (-1=detect, 0=lat/lon [deg], 1=UTM [m]).
Definition: mptrac.h:2637
ctl_t::qnt_idx
int qnt_idx
Quantity array index for air parcel IDs.
Definition: mptrac.h:2306
dd_t
Domain decomposition data structure.
Definition: mptrac.h:3772
met_t
Meteo data structure.
Definition: mptrac.h:3595
met_t::nx
int nx
Number of longitudes.
Definition: mptrac.h:3604
met_t::ny
int ny
Number of latitudes.
Definition: mptrac.h:3607
met_t::ps
float ps[EX][EY]
Surface pressure [hPa].
Definition: mptrac.h:3637
met_t::np
int np
Number of pressure levels.
Definition: mptrac.h:3610
met_t::pbl
float pbl[EX][EY]
Boundary layer pressure [hPa].
Definition: mptrac.h:3667
met_t::p
double p[EP]
Pressure levels [hPa].
Definition: mptrac.h:3622
Generated by doxygen 1.9.4