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 202 of file atm_init.c.

203 {
204
205 printf("\nMPTRAC atm_init tool.\n\n");
206 printf
207 ("Create an atmospheric data file with initial air parcel positions.\n");
208 printf("\n");
209 printf("Usage:\n");
210 printf(" atm_init <ctl> <atm_out> [KEY VALUE ...]\n");
211 printf("\n");
212 printf("Arguments:\n");
213 printf(" <ctl> Control file.\n");
214 printf(" <atm_out> Atmospheric output file.\n");
215 printf(" [KEY VALUE] Optional control parameters.\n");
216 printf("\nFurther information:\n");
217 printf(" Manual: https://slcs-jsc.github.io/mptrac/\n");
218}

◆ main()

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

Definition at line 39 of file atm_init.c.

41 {
42
43 atm_t *atm;
44
45 ctl_t ctl;
46
47 /* Allocate... */
48 ALLOC(atm, atm_t, 1);
49
50 /* Print usage information... */
51 USAGE;
52
53 /* Check arguments... */
54 if (argc < 3)
55 ERRMSG("Missing or invalid command-line arguments.\n\n"
56 "Usage: atm_init <ctl> <atm_out> [KEY VALUE ...]\n\n"
57 "Use -h for full help.");
58
59 /* Read control parameters... */
60 mptrac_read_ctl(argv[1], argc, argv, &ctl);
61 const int ens =
62 (int) scan_ctl(argv[1], argc, argv, "INIT_ENS", -1, "0", NULL);
63 const double t0 = scan_ctl(argv[1], argc, argv, "INIT_T0", -1, "0", NULL);
64 const double t1 = scan_ctl(argv[1], argc, argv, "INIT_T1", -1, "0", NULL);
65 const double dt = scan_ctl(argv[1], argc, argv, "INIT_DT", -1, "1", NULL);
66 const double z0 = scan_ctl(argv[1], argc, argv, "INIT_Z0", -1, "0", NULL);
67 const double z1 = scan_ctl(argv[1], argc, argv, "INIT_Z1", -1, "0", NULL);
68 const double dz = scan_ctl(argv[1], argc, argv, "INIT_DZ", -1, "1", NULL);
69 const double lon0 =
70 scan_ctl(argv[1], argc, argv, "INIT_LON0", -1, "0", NULL);
71 const double lon1 =
72 scan_ctl(argv[1], argc, argv, "INIT_LON1", -1, "0", NULL);
73 const double dlon =
74 scan_ctl(argv[1], argc, argv, "INIT_DLON", -1, "1", NULL);
75 const double lat0 =
76 scan_ctl(argv[1], argc, argv, "INIT_LAT0", -1, "0", NULL);
77 const double lat1 =
78 scan_ctl(argv[1], argc, argv, "INIT_LAT1", -1, "0", NULL);
79 const double dlat =
80 scan_ctl(argv[1], argc, argv, "INIT_DLAT", -1, "1", NULL);
81 const double st = scan_ctl(argv[1], argc, argv, "INIT_ST", -1, "0", NULL);
82 const double sz = scan_ctl(argv[1], argc, argv, "INIT_SZ", -1, "0", NULL);
83 const double slon =
84 scan_ctl(argv[1], argc, argv, "INIT_SLON", -1, "0", NULL);
85 const double slat =
86 scan_ctl(argv[1], argc, argv, "INIT_SLAT", -1, "0", NULL);
87 const double sx = scan_ctl(argv[1], argc, argv, "INIT_SX", -1, "0", NULL);
88 const double ut = scan_ctl(argv[1], argc, argv, "INIT_UT", -1, "0", NULL);
89 const double uz = scan_ctl(argv[1], argc, argv, "INIT_UZ", -1, "0", NULL);
90 const double ulon =
91 scan_ctl(argv[1], argc, argv, "INIT_ULON", -1, "0", NULL);
92 const double ulat =
93 scan_ctl(argv[1], argc, argv, "INIT_ULAT", -1, "0", NULL);
94 const int even =
95 (int) scan_ctl(argv[1], argc, argv, "INIT_EVENLY", -1, "0", NULL);
96 const int rep =
97 (int) scan_ctl(argv[1], argc, argv, "INIT_REP", -1, "1", NULL);
98 const double m = scan_ctl(argv[1], argc, argv, "INIT_MASS", -1, "0", NULL);
99 const double vmr = scan_ctl(argv[1], argc, argv, "INIT_VMR", -1, "0", NULL);
100 const double bellrad =
101 scan_ctl(argv[1], argc, argv, "INIT_BELLRAD", -1, "0", NULL);
102 const int idx_offset =
103 (int) scan_ctl(argv[1], argc, argv, "INIT_IDX_OFFSET", -1, "0", NULL);
104
105 /* Initialize random number generator... */
106 gsl_rng_env_setup();
107 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
108
109 /* Create grid... */
110 for (double t = t0; t <= t1; t += dt)
111 for (double z = z0; z <= z1; z += dz)
112 for (double lon = lon0; lon <= lon1; lon += dlon)
113 for (double lat = lat0; lat <= lat1; lat += dlat)
114 for (int irep = 0; irep < rep; irep++) {
115
116 /* Set position... */
117 double rg = gsl_ran_gaussian_ziggurat(rng, st / 2.3548);
118 double ru = ut * (gsl_rng_uniform(rng) - 0.5);
119 atm->time[atm->np] = (t + rg + ru);
120
121 rg = gsl_ran_gaussian_ziggurat(rng, sz / 2.3548);
122 ru = uz * (gsl_rng_uniform(rng) - 0.5);
123 atm->p[atm->np] = P(z + rg + ru);
124
125 rg = gsl_ran_gaussian_ziggurat(rng, slon / 2.3548);
126 double rx =
127 gsl_ran_gaussian_ziggurat(rng, DX2DEG(sx, lat) / 2.3548);
128 ru = ulon * (gsl_rng_uniform(rng) - 0.5);
129 atm->lon[atm->np] = (lon + rg + rx + ru);
130
131 if (ctl.qnt_ens >= 0)
132 atm->q[ctl.qnt_ens][atm->np] = ens;
133
134 do {
135 rg = gsl_ran_gaussian_ziggurat(rng, slat / 2.3548);
136 rx = gsl_ran_gaussian_ziggurat(rng, DY2DEG(sx) / 2.3548);
137 ru = ulat * (gsl_rng_uniform(rng) - 0.5);
138 atm->lat[atm->np] = (lat + rg + rx + ru);
139 } while (even && gsl_rng_uniform(rng) >
140 fabs(cos(DEG2RAD(atm->lat[atm->np]))));
141
142 /* Apply cosine bell (Williamson et al., 1992)... */
143 if (bellrad > 0) {
144 double x0[3], x1[3];
145 geo2cart(0.0, 0.5 * (lon0 + lon1), 0.5 * (lat0 + lat1), x0);
146 geo2cart(0.0, atm->lon[atm->np], atm->lat[atm->np], x1);
147 const double rad =
148 RE * acos(DOTP(x0, x1) / sqrt(DOTP(x0, x0)) /
149 sqrt(DOTP(x1, x1)));
150 if (rad > bellrad)
151 continue;
152 if (ctl.qnt_m >= 0)
153 atm->q[ctl.qnt_m][atm->np] =
154 0.5 * (1. + cos(M_PI * rad / bellrad));
155 if (ctl.qnt_vmr >= 0)
156 atm->q[ctl.qnt_vmr][atm->np] =
157 0.5 * (1. + cos(M_PI * rad / bellrad));
158 }
159
160 /* Set particle counter... */
161 if ((++atm->np) > NP)
162 ERRMSG("Too many particles!");
163 }
164
165 /* Check number of air parcels... */
166 if (atm->np <= 0)
167 ERRMSG("Did not create any air parcels!");
168
169 /* Initialize mass... */
170 if (ctl.qnt_m >= 0 && bellrad <= 0)
171 for (int ip = 0; ip < atm->np; ip++)
172 atm->q[ctl.qnt_m][ip] = m / atm->np;
173
174 /* Initialize volume mixing ratio... */
175 if (ctl.qnt_vmr >= 0 && bellrad <= 0)
176 for (int ip = 0; ip < atm->np; ip++)
177 atm->q[ctl.qnt_vmr][ip] = vmr;
178
179 /* Initialize air parcel index... */
180 if (ctl.qnt_idx >= 0)
181 for (int ip = 0; ip < atm->np; ip++)
182 atm->q[ctl.qnt_idx][ip] = idx_offset + ip;
183
184 /* Initialize age of air... */
185 if (ctl.qnt_aoa >= 0)
186 for (int ip = 0; ip < atm->np; ip++)
187 atm->q[ctl.qnt_aoa][ip] = atm->time[ip];
188
189 /* Save data... */
190 mptrac_write_atm(argv[2], &ctl, atm, 0);
191
192 /* Free... */
193 gsl_rng_free(rng);
194 free(atm);
195
196 return EXIT_SUCCESS;
197}
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:6526
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:10745
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:5248
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:1674
RE
#define RE
Mean radius of Earth [km].
Definition: mptrac.h:311
DOTP
#define DOTP(a, b)
Calculate the dot product of two vectors.
Definition: mptrac.h:736
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:2102
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: mptrac.h:1909
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1480
DX2DEG
#define DX2DEG(dx, lat)
Convert a distance in kilometers to degrees longitude at a given latitude.
Definition: mptrac.h:651
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:453
DEG2RAD
#define DEG2RAD(deg)
Converts degrees to radians.
Definition: mptrac.h:604
NP
#define NP
Maximum number of atmospheric data points.
Definition: mptrac.h:355
DY2DEG
#define DY2DEG(dy)
Convert a distance in kilometers to degrees latitude.
Definition: mptrac.h:669
atm_t
Air parcel data.
Definition: mptrac.h:3241
atm_t::time
double time[NP]
Time [s].
Definition: mptrac.h:3247
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3256
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3253
atm_t::np
int np
Number of air parcels.
Definition: mptrac.h:3244
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3259
atm_t::p
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3250
ctl_t
Control parameters.
Definition: mptrac.h:2190
ctl_t::qnt_m
int qnt_m
Quantity array index for mass.
Definition: mptrac.h:2221
ctl_t::qnt_aoa
int qnt_aoa
Quantity array index for age of air.
Definition: mptrac.h:2491
ctl_t::qnt_vmr
int qnt_vmr
Quantity array index for volume mixing ratio.
Definition: mptrac.h:2224
ctl_t::qnt_ens
int qnt_ens
Quantity array index for ensemble IDs.
Definition: mptrac.h:2215
ctl_t::qnt_idx
int qnt_idx
Quantity array index for air parcel IDs.
Definition: mptrac.h:2212
Here is the call graph for this function:
Generated by doxygen 1.9.4