41 {
42
44
46
48
50
51 FILE *out;
52
54
55 static float u[16], v[16], w[16];
56
58
59
64
65
67
68
69 if (argc < 4 && argc % 2 != 0)
70 ERRMSG(
"Missing or invalid command-line arguments.\n\n"
71 "Usage: met_subgrid <ctl> <subgrid.tab> <met0> <met1> [<met0b> <met1b> ...]\n\n"
72 "Use -h for full help.");
73
74
76
77
79
80
81 for (int i = 3; i < argc - 1; i += 2) {
82
83
85 ERRMSG(
"Cannot open file!");
87 ERRMSG(
"Cannot open file!");
88
89
90 for (
int ix = 0; ix < met0->
nx - 1; ix++)
91 for (
int iy = 0; iy < met0->
ny - 1; iy++)
92 for (
int iz = 0; iz < met0->
np - 1; iz++) {
93
94
95 u[0] = met0->
u[ix][iy][iz];
96 u[1] = met0->
u[ix + 1][iy][iz];
97 u[2] = met0->
u[ix][iy + 1][iz];
98 u[3] = met0->
u[ix + 1][iy + 1][iz];
99 u[4] = met0->
u[ix][iy][iz + 1];
100 u[5] = met0->
u[ix + 1][iy][iz + 1];
101 u[6] = met0->
u[ix][iy + 1][iz + 1];
102 u[7] = met0->
u[ix + 1][iy + 1][iz + 1];
103
104 v[0] = met0->
v[ix][iy][iz];
105 v[1] = met0->
v[ix + 1][iy][iz];
106 v[2] = met0->
v[ix][iy + 1][iz];
107 v[3] = met0->
v[ix + 1][iy + 1][iz];
108 v[4] = met0->
v[ix][iy][iz + 1];
109 v[5] = met0->
v[ix + 1][iy][iz + 1];
110 v[6] = met0->
v[ix][iy + 1][iz + 1];
111 v[7] = met0->
v[ix + 1][iy + 1][iz + 1];
112
113 w[0] = (float) (1e3 *
DP2DZ(met0->
w[ix][iy][iz], met0->
p[iz]));
114 w[1] = (float) (1e3 *
DP2DZ(met0->
w[ix + 1][iy][iz], met0->
p[iz]));
115 w[2] = (float) (1e3 *
DP2DZ(met0->
w[ix][iy + 1][iz], met0->
p[iz]));
116 w[3] =
117 (float) (1e3 *
DP2DZ(met0->
w[ix + 1][iy + 1][iz], met0->
p[iz]));
118 w[4] =
119 (float) (1e3 *
DP2DZ(met0->
w[ix][iy][iz + 1], met0->
p[iz + 1]));
120 w[5] =
121 (float) (1e3 *
122 DP2DZ(met0->
w[ix + 1][iy][iz + 1], met0->
p[iz + 1]));
123 w[6] =
124 (float) (1e3 *
125 DP2DZ(met0->
w[ix][iy + 1][iz + 1], met0->
p[iz + 1]));
126 w[7] =
127 (float) (1e3 *
128 DP2DZ(met0->
w[ix + 1][iy + 1][iz + 1], met0->
p[iz + 1]));
129
130
131 u[8] = met1->
u[ix][iy][iz];
132 u[9] = met1->
u[ix + 1][iy][iz];
133 u[10] = met1->
u[ix][iy + 1][iz];
134 u[11] = met1->
u[ix + 1][iy + 1][iz];
135 u[12] = met1->
u[ix][iy][iz + 1];
136 u[13] = met1->
u[ix + 1][iy][iz + 1];
137 u[14] = met1->
u[ix][iy + 1][iz + 1];
138 u[15] = met1->
u[ix + 1][iy + 1][iz + 1];
139
140 v[8] = met1->
v[ix][iy][iz];
141 v[9] = met1->
v[ix + 1][iy][iz];
142 v[10] = met1->
v[ix][iy + 1][iz];
143 v[11] = met1->
v[ix + 1][iy + 1][iz];
144 v[12] = met1->
v[ix][iy][iz + 1];
145 v[13] = met1->
v[ix + 1][iy][iz + 1];
146 v[14] = met1->
v[ix][iy + 1][iz + 1];
147 v[15] = met1->
v[ix + 1][iy + 1][iz + 1];
148
149 w[8] = (float) (1e3 *
DP2DZ(met1->
w[ix][iy][iz], met1->
p[iz]));
150 w[9] = (float) (1e3 *
DP2DZ(met1->
w[ix + 1][iy][iz], met1->
p[iz]));
151 w[10] = (float) (1e3 *
DP2DZ(met1->
w[ix][iy + 1][iz], met1->
p[iz]));
152 w[11] =
153 (float) (1e3 *
DP2DZ(met1->
w[ix + 1][iy + 1][iz], met1->
p[iz]));
154 w[12] =
155 (float) (1e3 *
DP2DZ(met1->
w[ix][iy][iz + 1], met1->
p[iz + 1]));
156 w[13] =
157 (float) (1e3 *
158 DP2DZ(met1->
w[ix + 1][iy][iz + 1], met1->
p[iz + 1]));
159 w[14] =
160 (float) (1e3 *
161 DP2DZ(met1->
w[ix][iy + 1][iz + 1], met1->
p[iz + 1]));
162 w[15] =
163 (float) (1e3 *
164 DP2DZ(met1->
w[ix + 1][iy + 1][iz + 1], met1->
p[iz + 1]));
165
166
167 usig[iz][iy] +=
stddev(u, 16);
168 vsig[iz][iy] +=
stddev(v, 16);
169 wsig[iz][iy] +=
stddev(w, 16);
170 n[iz][iy]++;
171
172
173 if (met0->
p[iz] > met0->
ps[ix][iy]
174 || met1->
p[iz] > met1->
ps[ix][iy]) {
175 usig[iz][iy] = NAN;
176 vsig[iz][iy] = NAN;
177 wsig[iz][iy] = NAN;
178 n[iz][iy] = 0;
179 }
180 }
181 }
182
183
184 LOG(1,
"Write subgrid data file: %s", argv[2]);
185 if (!(out = fopen(argv[2], "w")))
186 ERRMSG(
"Cannot create file!");
187
188
189 fprintf(out,
190 "# $1 = time [s]\n"
191 "# $2 = altitude [km]\n"
192 "# $3 = longitude [deg]\n"
193 "# $4 = latitude [deg]\n"
194 "# $5 = zonal wind standard deviation [m/s]\n"
195 "# $6 = meridional wind standard deviation [m/s]\n"
196 "# $7 = vertical velocity standard deviation [m/s]\n"
197 "# $8 = number of data points\n");
198
199
200 for (
int iy = 0; iy < met0->
ny - 1; iy++) {
201 fprintf(out, "\n");
202 for (
int iz = 0; iz < met0->
np - 1; iz++)
203 fprintf(out, "%.2f %g %g %g %g %g %g %d\n",
205 0.5 * (
Z(met0->
p[iz]) +
Z(met1->
p[iz + 1])),
206 0.0, 0.5 * (met0->
lat[iy] + met1->
lat[iy + 1]),
207 usig[iz][iy] / n[iz][iy], vsig[iz][iy] / n[iz][iy],
208 wsig[iz][iy] / n[iz][iy], n[iz][iy]);
209 }
210
211
212 fclose(out);
213
214
215 free(clim);
216 free(met0);
217 free(met1);
218 free(dd);
219
220 return EXIT_SUCCESS;
221}
float stddev(const float *data, const int n)
Calculates the standard deviation of a set of data.
void mptrac_read_clim(const ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
int mptrac_read_met(const char *filename, const ctl_t *ctl, const clim_t *clim, met_t *met, dd_t *dd)
Reads meteorological data from a file, supporting multiple formats and MPI broadcasting.
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.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
#define EY
Maximum number of latitudes for meteo data.
#define USAGE
Print usage information on -h or --help.
#define Z(p)
Convert pressure to altitude.
#define DP2DZ(dp, p)
Convert a pressure difference to a height difference in the vertical direction.
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
#define LOG(level,...)
Print a log message with a specified logging level.
#define EP
Maximum number of pressure levels for meteo data.
Domain decomposition data structure.
float w[EX][EY][EP]
Vertical velocity [hPa/s].
int nx
Number of longitudes.
int ny
Number of latitudes.
float ps[EX][EY]
Surface pressure [hPa].
int np
Number of pressure levels.
float u[EX][EY][EP]
Zonal wind [m/s].
float v[EX][EY][EP]
Meridional wind [m/s].
double lat[EY]
Latitudes [deg].
double p[EP]
Pressure levels [hPa].