MPTRAC
Functions
atm_init.c File Reference

Create atmospheric data file with initial air parcel positions. More...

#include "mptrac.h"

Go to the source code of this file.

Functions

void usage (void)
 Print command-line help. More...
 
int main (int argc, char *argv[])
 

Detailed Description

Create atmospheric data file with initial air parcel positions.

Definition in file atm_init.c.

Function Documentation

◆ usage()

void usage ( void  )

Print command-line help.

Definition at line 259 of file atm_init.c.

260 {
261
262 printf("\nMPTRAC atm_init tool.\n\n");
263 printf
264 ("Create an atmospheric data file with initial air parcel positions.\n");
265 printf("\n");
266 printf("Usage:\n");
267 printf(" atm_init <ctl> <atm_out> [KEY VALUE ...]\n");
268 printf("\n");
269 printf("Arguments:\n");
270 printf(" <ctl> Control file.\n");
271 printf(" <atm_out> Atmospheric output file.\n");
272 printf(" [KEY VALUE] Optional control parameters.\n");
273 printf("\nFurther information:\n");
274 printf(" Manual: https://slcs-jsc.github.io/mptrac/\n");
275}

◆ main()

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

Definition at line 39 of file atm_init.c.

41 {
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 lon0 =
80 scan_ctl(argv[1], argc, argv, "INIT_LON0", -1, "0", NULL);
81 const double lon1 =
82 scan_ctl(argv[1], argc, argv, "INIT_LON1", -1, "0", NULL);
83 const double dlon =
84 scan_ctl(argv[1], argc, argv, "INIT_DLON", -1, "1", NULL);
85 const double lat0 =
86 scan_ctl(argv[1], argc, argv, "INIT_LAT0", -1, "0", NULL);
87 const double lat1 =
88 scan_ctl(argv[1], argc, argv, "INIT_LAT1", -1, "0", NULL);
89 const double dlat =
90 scan_ctl(argv[1], argc, argv, "INIT_DLAT", -1, "1", NULL);
91 const double st = scan_ctl(argv[1], argc, argv, "INIT_ST", -1, "0", NULL);
92 const double sz = scan_ctl(argv[1], argc, argv, "INIT_SZ", -1, "0", NULL);
93 const double slon =
94 scan_ctl(argv[1], argc, argv, "INIT_SLON", -1, "0", NULL);
95 const double slat =
96 scan_ctl(argv[1], argc, argv, "INIT_SLAT", -1, "0", NULL);
97 const double sx = scan_ctl(argv[1], argc, argv, "INIT_SX", -1, "0", NULL);
98 const double ut = scan_ctl(argv[1], argc, argv, "INIT_UT", -1, "0", NULL);
99 const double uz = scan_ctl(argv[1], argc, argv, "INIT_UZ", -1, "0", NULL);
100 const double ulon =
101 scan_ctl(argv[1], argc, argv, "INIT_ULON", -1, "0", NULL);
102 const double ulat =
103 scan_ctl(argv[1], argc, argv, "INIT_ULAT", -1, "0", NULL);
104 const int even =
105 (int) scan_ctl(argv[1], argc, argv, "INIT_EVENLY", -1, "0", NULL);
106 const int rho =
107 (int) scan_ctl(argv[1], argc, argv, "INIT_RHO", -1, "0", NULL);
108 const int rep =
109 (int) scan_ctl(argv[1], argc, argv, "INIT_REP", -1, "1", NULL);
110 const double m = scan_ctl(argv[1], argc, argv, "INIT_MASS", -1, "0", NULL);
111 const double vmr = scan_ctl(argv[1], argc, argv, "INIT_VMR", -1, "0", NULL);
112 const double bellrad =
113 scan_ctl(argv[1], argc, argv, "INIT_BELLRAD", -1, "0", NULL);
114 const int idx_offset =
115 (int) scan_ctl(argv[1], argc, argv, "INIT_IDX_OFFSET", -1, "0", NULL);
116
117 /* Check arguments... */
118 if (ctl.met_coord_type != 0 && even)
119 ERRMSG("INIT_EVENLY is only supported for lat/lon grids");
120 if (ctl.met_coord_type != 0 && bellrad > 0)
121 ERRMSG("INIT_BELLRAD is only supported for lat/lon grids");
122
123 /* Read climatological data... */
124 mptrac_read_clim(&ctl, clim);
125
126 /* Initialize random number generator... */
127 gsl_rng_env_setup();
128 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
129
130 /* Get meteorological data for density-weighted initialization... */
131 double ptop, ztop, rhomax = RHO(1100., 200.);
132 if (rho) {
133 mptrac_get_met(&ctl, clim, t0, &met0, &met1, dd);
134 ptop = gsl_stats_min(met0->p, 1, (size_t)met0->np);
135 ztop = Z(ptop);
136 }
137
138 /* Create grid... */
139 for (double t = t0; t <= t1; t += dt)
140 for (double z = z0; z <= z1; z += dz)
141 for (double lon = lon0; lon <= lon1; lon += dlon)
142 for (double lat = lat0; lat <= lat1; lat += dlat)
143 for (int irep = 0; irep < rep; irep++) {
144
145 /* Set time... */
146 double rg = gsl_ran_gaussian_ziggurat(rng, st / 2.3548);
147 double ru = ut * (gsl_rng_uniform(rng) - 0.5);
148 atm->time[atm->np] = (t + rg + ru);
149
150 /* Set pressure... */
151 rg = gsl_ran_gaussian_ziggurat(rng, sz / 2.3548);
152 ru = uz * (gsl_rng_uniform(rng) - 0.5);
153 atm->p[atm->np] = P(z + rg + ru);
154
155 /* Set horizontal position... */
156 rg = gsl_ran_gaussian_ziggurat(rng, slon / 2.3548);
157 const double sx_coord =
158 ctl.met_coord_type == 0 ? DX2DEG(sx, lat) : sx;
159 double rx = gsl_ran_gaussian_ziggurat(rng, sx_coord / 2.3548);
160 ru = ulon * (gsl_rng_uniform(rng) - 0.5);
161 atm->lon[atm->np] = (lon + rg + rx + ru);
162
163 /* Set ensemble index... */
164 if (ctl.qnt_ens >= 0)
165 atm->q[ctl.qnt_ens][atm->np] = ens;
166
167 /* Apply cosine-latitude weighting... */
168 do {
169 rg = gsl_ran_gaussian_ziggurat(rng, slat / 2.3548);
170 const double sy_coord =
171 ctl.met_coord_type == 0 ? DY2DEG(sx) : sx;
172 rx = gsl_ran_gaussian_ziggurat(rng, sy_coord / 2.3548);
173 ru = ulat * (gsl_rng_uniform(rng) - 0.5);
174 atm->lat[atm->np] = (lat + rg + rx + ru);
175 } while (even && gsl_rng_uniform(rng) >
176 fabs(cos(DEG2RAD(atm->lat[atm->np]))));
177
178 /* Apply cosine bell (Williamson et al., 1992)... */
179 if (bellrad > 0) {
180 double x0[3], x1[3];
181 geo2cart(0.0, 0.5 * (lon0 + lon1), 0.5 * (lat0 + lat1), x0);
182 geo2cart(0.0, atm->lon[atm->np], atm->lat[atm->np], x1);
183 const double rad =
184 RE * acos(DOTP(x0, x1) / sqrt(DOTP(x0, x0)) /
185 sqrt(DOTP(x1, x1)));
186 if (rad > bellrad)
187 continue;
188 if (ctl.qnt_m >= 0)
189 atm->q[ctl.qnt_m][atm->np] =
190 0.5 * (1. + cos(M_PI * rad / bellrad));
191 if (ctl.qnt_vmr >= 0)
192 atm->q[ctl.qnt_vmr][atm->np] =
193 0.5 * (1. + cos(M_PI * rad / bellrad));
194 }
195
196 /* Apply density weighting... */
197 if (rho) {
198 INTPOL_INIT;
199 double ps, ttry;
200 intpol_met_time_2d(met0, met0->ps, met1, met1->ps,
201 atm->time[atm->np], atm->lon[atm->np],
202 atm->lat[atm->np], &ps, ci, cw, 1);
203 const double zs = Z(ps);
204 do {
205 atm->p[atm->np] = P(zs + (ztop - zs) * gsl_rng_uniform(rng));
206 intpol_met_time_3d(met0, met0->t, met1, met1->t,
207 atm->time[atm->np], atm->p[atm->np],
208 atm->lon[atm->np], atm->lat[atm->np],
209 &ttry, ci, cw, 0);
210 } while (RHO(atm->p[atm->np], ttry) <= rhomax);
211 }
212
213 /* Set particle counter... */
214 if ((++atm->np) > NP)
215 ERRMSG("Too many particles!");
216 }
217
218 /* Check number of air parcels... */
219 if (atm->np <= 0)
220 ERRMSG("Did not create any air parcels!");
221
222 /* Initialize mass... */
223 if (ctl.qnt_m >= 0 && bellrad <= 0)
224 for (int ip = 0; ip < atm->np; ip++)
225 atm->q[ctl.qnt_m][ip] = m / atm->np;
226
227 /* Initialize volume mixing ratio... */
228 if (ctl.qnt_vmr >= 0 && bellrad <= 0)
229 for (int ip = 0; ip < atm->np; ip++)
230 atm->q[ctl.qnt_vmr][ip] = vmr;
231
232 /* Initialize air parcel index... */
233 if (ctl.qnt_idx >= 0)
234 for (int ip = 0; ip < atm->np; ip++)
235 atm->q[ctl.qnt_idx][ip] = idx_offset + ip;
236
237 /* Initialize age of air... */
238 if (ctl.qnt_aoa >= 0)
239 for (int ip = 0; ip < atm->np; ip++)
240 atm->q[ctl.qnt_aoa][ip] = atm->time[ip];
241
242 /* Save data... */
243 mptrac_write_atm(argv[2], &ctl, atm, 0);
244
245 /* Free... */
246 gsl_rng_free(rng);
247 free(clim);
248 free(atm);
249 free(met0);
250 free(met1);
251 free(dd);
252
253 return EXIT_SUCCESS;
254}
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:7553
intpol_met_time_3d
void intpol_met_time_3d(const met_t *met0, float array0[EX][EY][EP], const met_t *met1, float array1[EX][EY][EP], const double ts, const double p, const double lon, const double lat, double *var, int *ci, double *cw, const int init)
Interpolates meteorological data in 3D space and time.
Definition: mptrac.c:3059
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:11856
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:5944
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:6163
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:3088
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:6223
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:2548
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:802
INTPOL_INIT
#define INTPOL_INIT
Initialize arrays for interpolation.
Definition: mptrac.h:940
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:2172
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: mptrac.h:1979
Z
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:2009
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1550
DX2DEG
#define DX2DEG(dx, lat)
Convert a distance in kilometers to degrees longitude at a given latitude.
Definition: mptrac.h:670
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:623
NP
#define NP
Maximum number of atmospheric data points.
Definition: mptrac.h:359
RHO
#define RHO(p, t)
Compute density of air.
Definition: mptrac.h:1727
DY2DEG
#define DY2DEG(dy)
Convert a distance in kilometers to degrees latitude.
Definition: mptrac.h:688
atm_t
Air parcel data.
Definition: mptrac.h:3311
atm_t::time
double time[NP]
Time [s].
Definition: mptrac.h:3317
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3326
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3323
atm_t::np
int np
Number of air parcels.
Definition: mptrac.h:3314
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3329
atm_t::p
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3320
clim_t
Climatological data.
Definition: mptrac.h:3506
ctl_t
Control parameters.
Definition: mptrac.h:2260
ctl_t::qnt_m
int qnt_m
Quantity array index for mass.
Definition: mptrac.h:2291
ctl_t::qnt_aoa
int qnt_aoa
Quantity array index for age of air.
Definition: mptrac.h:2561
ctl_t::qnt_vmr
int qnt_vmr
Quantity array index for volume mixing ratio.
Definition: mptrac.h:2294
ctl_t::qnt_ens
int qnt_ens
Quantity array index for ensemble IDs.
Definition: mptrac.h:2285
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:2613
ctl_t::qnt_idx
int qnt_idx
Quantity array index for air parcel IDs.
Definition: mptrac.h:2282
dd_t
Domain decomposition data structure.
Definition: mptrac.h:3742
met_t
Meteo data structure.
Definition: mptrac.h:3565
met_t::ps
float ps[EX][EY]
Surface pressure [hPa].
Definition: mptrac.h:3607
met_t::np
int np
Number of pressure levels.
Definition: mptrac.h:3580
met_t::t
float t[EX][EY][EP]
Temperature [K].
Definition: mptrac.h:3682
met_t::p
double p[EP]
Pressure levels [hPa].
Definition: mptrac.h:3592
Here is the call graph for this function:
Generated by doxygen 1.9.4