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

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

◆ main()

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

Definition at line 27 of file atm_init.c.

29 {
30
31 atm_t *atm;
32
33 ctl_t ctl;
34
35 gsl_rng *rng;
36
37 /* Allocate... */
38 ALLOC(atm, atm_t, 1);
39
40 /* Check arguments... */
41 if (argc < 3)
42 ERRMSG("Give parameters: <ctl> <atm_out>");
43
44 /* Read control parameters... */
45 read_ctl(argv[1], argc, argv, &ctl);
46 double t0 = scan_ctl(argv[1], argc, argv, "INIT_T0", -1, "0", NULL);
47 double t1 = scan_ctl(argv[1], argc, argv, "INIT_T1", -1, "0", NULL);
48 double dt = scan_ctl(argv[1], argc, argv, "INIT_DT", -1, "1", NULL);
49 double z0 = scan_ctl(argv[1], argc, argv, "INIT_Z0", -1, "0", NULL);
50 double z1 = scan_ctl(argv[1], argc, argv, "INIT_Z1", -1, "0", NULL);
51 double dz = scan_ctl(argv[1], argc, argv, "INIT_DZ", -1, "1", NULL);
52 double lon0 = scan_ctl(argv[1], argc, argv, "INIT_LON0", -1, "0", NULL);
53 double lon1 = scan_ctl(argv[1], argc, argv, "INIT_LON1", -1, "0", NULL);
54 double dlon = scan_ctl(argv[1], argc, argv, "INIT_DLON", -1, "1", NULL);
55 double lat0 = scan_ctl(argv[1], argc, argv, "INIT_LAT0", -1, "0", NULL);
56 double lat1 = scan_ctl(argv[1], argc, argv, "INIT_LAT1", -1, "0", NULL);
57 double dlat = scan_ctl(argv[1], argc, argv, "INIT_DLAT", -1, "1", NULL);
58 double st = scan_ctl(argv[1], argc, argv, "INIT_ST", -1, "0", NULL);
59 double sz = scan_ctl(argv[1], argc, argv, "INIT_SZ", -1, "0", NULL);
60 double slon = scan_ctl(argv[1], argc, argv, "INIT_SLON", -1, "0", NULL);
61 double slat = scan_ctl(argv[1], argc, argv, "INIT_SLAT", -1, "0", NULL);
62 double sx = scan_ctl(argv[1], argc, argv, "INIT_SX", -1, "0", NULL);
63 double ut = scan_ctl(argv[1], argc, argv, "INIT_UT", -1, "0", NULL);
64 double uz = scan_ctl(argv[1], argc, argv, "INIT_UZ", -1, "0", NULL);
65 double ulon = scan_ctl(argv[1], argc, argv, "INIT_ULON", -1, "0", NULL);
66 double ulat = scan_ctl(argv[1], argc, argv, "INIT_ULAT", -1, "0", NULL);
67 int even =
68 (int) scan_ctl(argv[1], argc, argv, "INIT_EVENLY", -1, "0", NULL);
69 int rep = (int) scan_ctl(argv[1], argc, argv, "INIT_REP", -1, "1", NULL);
70 double m = scan_ctl(argv[1], argc, argv, "INIT_MASS", -1, "0", NULL);
71 double vmr = scan_ctl(argv[1], argc, argv, "INIT_VMR", -1, "0", NULL);
72 double bellrad =
73 scan_ctl(argv[1], argc, argv, "INIT_BELLRAD", -1, "0", NULL);
74
75 /* Initialize random number generator... */
76 gsl_rng_env_setup();
77 rng = gsl_rng_alloc(gsl_rng_default);
78
79 /* Create grid... */
80 for (double t = t0; t <= t1; t += dt)
81 for (double z = z0; z <= z1; z += dz)
82 for (double lon = lon0; lon <= lon1; lon += dlon)
83 for (double lat = lat0; lat <= lat1; lat += dlat)
84 for (int irep = 0; irep < rep; irep++) {
85
86 /* Set position... */
87 double rg = gsl_ran_gaussian_ziggurat(rng, st / 2.3548);
88 double ru = ut * (gsl_rng_uniform(rng) - 0.5);
89 atm->time[atm->np] = (t + rg + ru);
90
91 rg = gsl_ran_gaussian_ziggurat(rng, sz / 2.3548);
92 ru = uz * (gsl_rng_uniform(rng) - 0.5);
93 atm->p[atm->np] = P(z + rg + ru);
94
95 rg = gsl_ran_gaussian_ziggurat(rng, slon / 2.3548);
96 double rx =
97 gsl_ran_gaussian_ziggurat(rng, DX2DEG(sx, lat) / 2.3548);
98 ru = ulon * (gsl_rng_uniform(rng) - 0.5);
99 atm->lon[atm->np] = (lon + rg + rx + ru);
100
101 do {
102 rg = gsl_ran_gaussian_ziggurat(rng, slat / 2.3548);
103 rx = gsl_ran_gaussian_ziggurat(rng, DY2DEG(sx) / 2.3548);
104 ru = ulat * (gsl_rng_uniform(rng) - 0.5);
105 atm->lat[atm->np] = (lat + rg + rx + ru);
106 } while (even && gsl_rng_uniform(rng) >
107 fabs(cos(atm->lat[atm->np] * M_PI / 180.)));
108
109 /* Apply cosine bell (Williamson et al., 1992)... */
110 if (bellrad > 0) {
111 double x0[3], x1[3];
112 geo2cart(0.0, 0.5 * (lon0 + lon1), 0.5 * (lat0 + lat1), x0);
113 geo2cart(0.0, atm->lon[atm->np], atm->lat[atm->np], x1);
114 double rad = RE * acos(DOTP(x0, x1) / NORM(x0) / NORM(x1));
115 if (rad > bellrad)
116 continue;
117 if (ctl.qnt_m >= 0)
118 atm->q[ctl.qnt_m][atm->np] =
119 0.5 * (1. + cos(M_PI * rad / bellrad));
120 if (ctl.qnt_vmr >= 0)
121 atm->q[ctl.qnt_vmr][atm->np] =
122 0.5 * (1. + cos(M_PI * rad / bellrad));
123 }
124
125 /* Set particle counter... */
126 if ((++atm->np) > NP)
127 ERRMSG("Too many particles!");
128 }
129
130 /* Check number of air parcels... */
131 if (atm->np <= 0)
132 ERRMSG("Did not create any air parcels!");
133
134 /* Initialize mass... */
135 if (ctl.qnt_m >= 0 && bellrad <= 0)
136 for (int ip = 0; ip < atm->np; ip++)
137 atm->q[ctl.qnt_m][ip] = m / atm->np;
138
139 /* Initialize volume mixing ratio... */
140 if (ctl.qnt_vmr >= 0 && bellrad <= 0)
141 for (int ip = 0; ip < atm->np; ip++)
142 atm->q[ctl.qnt_vmr][ip] = vmr;
143
144 /* Initialize air parcel index... */
145 if (ctl.qnt_idx >= 0)
146 for (int ip = 0; ip < atm->np; ip++)
147 atm->q[ctl.qnt_idx][ip] = ip;
148
149 /* Initialize age of air... */
150 if (ctl.qnt_aoa >= 0)
151 for (int ip = 0; ip < atm->np; ip++)
152 atm->q[ctl.qnt_aoa][ip] = atm->time[ip];
153
154 /* Save data... */
155 write_atm(argv[2], &ctl, atm, 0);
156
157 /* Free... */
158 gsl_rng_free(rng);
159 free(atm);
160
161 return EXIT_SUCCESS;
162}
write_atm
void write_atm(const char *filename, ctl_t *ctl, atm_t *atm, double t)
Writes air parcel data to a file in various formats.
Definition: mptrac.c:8317
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:4789
scan_ctl
double scan_ctl(const char *filename, int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Scans a control file or command-line arguments for a specified variable.
Definition: mptrac.c:7943
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:941
RE
#define RE
Mean radius of Earth [km].
Definition: mptrac.h:217
DOTP
#define DOTP(a, b)
Calculate the dot product of two vectors.
Definition: mptrac.h:584
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:1881
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1255
DX2DEG
#define DX2DEG(dx, lat)
Convert a distance in kilometers to degrees longitude at a given latitude.
Definition: mptrac.h:501
NORM
#define NORM(a)
Compute the norm (magnitude) of a vector.
Definition: mptrac.h:1198
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:344
NP
#define NP
Maximum number of atmospheric data points.
Definition: mptrac.h:241
DY2DEG
#define DY2DEG(dy)
Convert a distance in kilometers to degrees latitude.
Definition: mptrac.h:519
atm_t
Air parcel data.
Definition: mptrac.h:3068
atm_t::time
double time[NP]
Time [s].
Definition: mptrac.h:3074
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3083
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3080
atm_t::np
int np
Number of air parcels.
Definition: mptrac.h:3071
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3086
atm_t::p
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3077
ctl_t
Control parameters.
Definition: mptrac.h:2135
ctl_t::qnt_m
int qnt_m
Quantity array index for mass.
Definition: mptrac.h:2184
ctl_t::qnt_aoa
int qnt_aoa
Quantity array index for age of air.
Definition: mptrac.h:2436
ctl_t::qnt_vmr
int qnt_vmr
Quantity array index for volume mixing ratio.
Definition: mptrac.h:2187
ctl_t::qnt_idx
int qnt_idx
Quantity array index for air parcel IDs.
Definition: mptrac.h:2175
Generated on Mon Sep 16 2024 17:36:25 for MPTRAC by doxygen 1.9.4