29 {
30
32
34
35
37
38
39 if (argc < 3)
40 ERRMSG(
"Give parameters: <ctl> <atm_out>");
41
42
44 const int ens =
45 (int)
scan_ctl(argv[1], argc, argv,
"INIT_ENS", -1,
"0", NULL);
46 const double t0 =
scan_ctl(argv[1], argc, argv,
"INIT_T0", -1,
"0", NULL);
47 const double t1 =
scan_ctl(argv[1], argc, argv,
"INIT_T1", -1,
"0", NULL);
48 const double dt =
scan_ctl(argv[1], argc, argv,
"INIT_DT", -1,
"1", NULL);
49 const double z0 =
scan_ctl(argv[1], argc, argv,
"INIT_Z0", -1,
"0", NULL);
50 const double z1 =
scan_ctl(argv[1], argc, argv,
"INIT_Z1", -1,
"0", NULL);
51 const double dz =
scan_ctl(argv[1], argc, argv,
"INIT_DZ", -1,
"1", NULL);
52 const double lon0 =
53 scan_ctl(argv[1], argc, argv,
"INIT_LON0", -1,
"0", NULL);
54 const double lon1 =
55 scan_ctl(argv[1], argc, argv,
"INIT_LON1", -1,
"0", NULL);
56 const double dlon =
57 scan_ctl(argv[1], argc, argv,
"INIT_DLON", -1,
"1", NULL);
58 const double lat0 =
59 scan_ctl(argv[1], argc, argv,
"INIT_LAT0", -1,
"0", NULL);
60 const double lat1 =
61 scan_ctl(argv[1], argc, argv,
"INIT_LAT1", -1,
"0", NULL);
62 const double dlat =
63 scan_ctl(argv[1], argc, argv,
"INIT_DLAT", -1,
"1", NULL);
64 const double st =
scan_ctl(argv[1], argc, argv,
"INIT_ST", -1,
"0", NULL);
65 const double sz =
scan_ctl(argv[1], argc, argv,
"INIT_SZ", -1,
"0", NULL);
66 const double slon =
67 scan_ctl(argv[1], argc, argv,
"INIT_SLON", -1,
"0", NULL);
68 const double slat =
69 scan_ctl(argv[1], argc, argv,
"INIT_SLAT", -1,
"0", NULL);
70 const double sx =
scan_ctl(argv[1], argc, argv,
"INIT_SX", -1,
"0", NULL);
71 const double ut =
scan_ctl(argv[1], argc, argv,
"INIT_UT", -1,
"0", NULL);
72 const double uz =
scan_ctl(argv[1], argc, argv,
"INIT_UZ", -1,
"0", NULL);
73 const double ulon =
74 scan_ctl(argv[1], argc, argv,
"INIT_ULON", -1,
"0", NULL);
75 const double ulat =
76 scan_ctl(argv[1], argc, argv,
"INIT_ULAT", -1,
"0", NULL);
77 const int even =
78 (int)
scan_ctl(argv[1], argc, argv,
"INIT_EVENLY", -1,
"0", NULL);
79 const int rep =
80 (int)
scan_ctl(argv[1], argc, argv,
"INIT_REP", -1,
"1", NULL);
81 const double m =
scan_ctl(argv[1], argc, argv,
"INIT_MASS", -1,
"0", NULL);
82 const double vmr =
scan_ctl(argv[1], argc, argv,
"INIT_VMR", -1,
"0", NULL);
83 const double bellrad =
84 scan_ctl(argv[1], argc, argv,
"INIT_BELLRAD", -1,
"0", NULL);
85
86
87 gsl_rng_env_setup();
88 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
89
90
91 for (double t = t0; t <= t1; t += dt)
92 for (double z = z0; z <= z1; z += dz)
93 for (double lon = lon0; lon <= lon1; lon += dlon)
94 for (double lat = lat0; lat <= lat1; lat += dlat)
95 for (int irep = 0; irep < rep; irep++) {
96
97
98 double rg = gsl_ran_gaussian_ziggurat(rng, st / 2.3548);
99 double ru = ut * (gsl_rng_uniform(rng) - 0.5);
100 atm->
time[atm->
np] = (t + rg + ru);
101
102 rg = gsl_ran_gaussian_ziggurat(rng, sz / 2.3548);
103 ru = uz * (gsl_rng_uniform(rng) - 0.5);
104 atm->
p[atm->
np] =
P(z + rg + ru);
105
106 rg = gsl_ran_gaussian_ziggurat(rng, slon / 2.3548);
107 double rx =
108 gsl_ran_gaussian_ziggurat(rng,
DX2DEG(sx, lat) / 2.3548);
109 ru = ulon * (gsl_rng_uniform(rng) - 0.5);
110 atm->
lon[atm->
np] = (lon + rg + rx + ru);
111
114
115 do {
116 rg = gsl_ran_gaussian_ziggurat(rng, slat / 2.3548);
117 rx = gsl_ran_gaussian_ziggurat(rng,
DY2DEG(sx) / 2.3548);
118 ru = ulat * (gsl_rng_uniform(rng) - 0.5);
119 atm->
lat[atm->
np] = (lat + rg + rx + ru);
120 } while (even && gsl_rng_uniform(rng) >
122
123
124 if (bellrad > 0) {
125 double x0[3], x1[3];
126 geo2cart(0.0, 0.5 * (lon0 + lon1), 0.5 * (lat0 + lat1), x0);
128 const double rad =
130 if (rad > bellrad)
131 continue;
134 0.5 * (1. + cos(M_PI * rad / bellrad));
137 0.5 * (1. + cos(M_PI * rad / bellrad));
138 }
139
140
141 if ((++atm->
np) >
NP)
142 ERRMSG(
"Too many particles!");
143 }
144
145
147 ERRMSG(
"Did not create any air parcels!");
148
149
150 if (ctl.
qnt_m >= 0 && bellrad <= 0)
151 for (
int ip = 0; ip < atm->
np; ip++)
152 atm->
q[ctl.
qnt_m][ip] = m / atm->
np;
153
154
155 if (ctl.
qnt_vmr >= 0 && bellrad <= 0)
156 for (
int ip = 0; ip < atm->
np; ip++)
158
159
161 for (
int ip = 0; ip < atm->
np; ip++)
163
164
166 for (
int ip = 0; ip < atm->
np; ip++)
168
169
171
172
173 gsl_rng_free(rng);
174 free(atm);
175
176 return EXIT_SUCCESS;
177}
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.
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_ctl(const char *filename, int argc, char *argv[], ctl_t *ctl)
Reads control parameters from a configuration file and populates the given structure.
void geo2cart(const double z, const double lon, const double lat, double *x)
Converts geographic coordinates (longitude, latitude, altitude) to Cartesian coordinates.
#define RE
Mean radius of Earth [km].
#define DOTP(a, b)
Calculate the dot product of two vectors.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
#define P(z)
Compute pressure at given altitude.
#define DX2DEG(dx, lat)
Convert a distance in kilometers to degrees longitude at a given latitude.
#define NORM(a)
Compute the norm (magnitude) of a vector.
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
#define DEG2RAD(deg)
Converts degrees to radians.
#define NP
Maximum number of atmospheric data points.
#define DY2DEG(dy)
Convert a distance in kilometers to degrees latitude.
double lat[NP]
Latitude [deg].
double lon[NP]
Longitude [deg].
int np
Number of air parcels.
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
double p[NP]
Pressure [hPa].
int qnt_m
Quantity array index for mass.
int qnt_aoa
Quantity array index for age of air.
int qnt_vmr
Quantity array index for volume mixing ratio.
int qnt_ens
Quantity array index for ensemble IDs.
int qnt_idx
Quantity array index for air parcel IDs.