MPTRAC
met_subgrid.c
Go to the documentation of this file.
1/*
2 This file is part of MPTRAC.
3
4 MPTRAC is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 MPTRAC is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
16
17 Copyright (C) 2013-2025 Forschungszentrum Juelich GmbH
18*/
19
25#include "mptrac.h"
26
27/* ------------------------------------------------------------
28 Main...
29 ------------------------------------------------------------ */
30
31int main(
32 int argc,
33 char *argv[]) {
34
35 ctl_t ctl;
36
37 clim_t *clim;
38
39 met_t *met0, *met1;
40
41 dd_t *dd;
42
43 FILE *out;
44
45 static double usig[EP][EY], vsig[EP][EY], wsig[EP][EY];
46
47 static float u[16], v[16], w[16];
48
49 static int n[EP][EY];
50
51 /* Allocate... */
52 ALLOC(clim, clim_t, 1);
53 ALLOC(met0, met_t, 1);
54 ALLOC(met1, met_t, 1);
55 ALLOC(dd, dd_t, 1);
56
57 /* Check arguments... */
58 if (argc < 4 && argc % 2 != 0)
59 ERRMSG
60 ("Give parameters: <ctl> <subgrid.tab> <met0> <met1> [ <met0> <met1> ... ]");
61
62 /* Read control parameters... */
63 mptrac_read_ctl(argv[1], argc, argv, &ctl);
64
65 /* Read climatological data... */
66 mptrac_read_clim(&ctl, clim);
67
68 /* Loop over data files... */
69 for (int i = 3; i < argc - 1; i += 2) {
70
71 /* Read meteorological data... */
72 if (!mptrac_read_met(argv[i], &ctl, clim, met0, dd))
73 ERRMSG("Cannot open file!");
74 if (!mptrac_read_met(argv[i + 1], &ctl, clim, met1, dd))
75 ERRMSG("Cannot open file!");
76
77 /* Loop over grid boxes... */
78 for (int ix = 0; ix < met0->nx - 1; ix++)
79 for (int iy = 0; iy < met0->ny - 1; iy++)
80 for (int iz = 0; iz < met0->np - 1; iz++) {
81
82 /* Collect local wind data... */
83 u[0] = met0->u[ix][iy][iz];
84 u[1] = met0->u[ix + 1][iy][iz];
85 u[2] = met0->u[ix][iy + 1][iz];
86 u[3] = met0->u[ix + 1][iy + 1][iz];
87 u[4] = met0->u[ix][iy][iz + 1];
88 u[5] = met0->u[ix + 1][iy][iz + 1];
89 u[6] = met0->u[ix][iy + 1][iz + 1];
90 u[7] = met0->u[ix + 1][iy + 1][iz + 1];
91
92 v[0] = met0->v[ix][iy][iz];
93 v[1] = met0->v[ix + 1][iy][iz];
94 v[2] = met0->v[ix][iy + 1][iz];
95 v[3] = met0->v[ix + 1][iy + 1][iz];
96 v[4] = met0->v[ix][iy][iz + 1];
97 v[5] = met0->v[ix + 1][iy][iz + 1];
98 v[6] = met0->v[ix][iy + 1][iz + 1];
99 v[7] = met0->v[ix + 1][iy + 1][iz + 1];
100
101 w[0] = (float) (1e3 * DP2DZ(met0->w[ix][iy][iz], met0->p[iz]));
102 w[1] = (float) (1e3 * DP2DZ(met0->w[ix + 1][iy][iz], met0->p[iz]));
103 w[2] = (float) (1e3 * DP2DZ(met0->w[ix][iy + 1][iz], met0->p[iz]));
104 w[3] =
105 (float) (1e3 * DP2DZ(met0->w[ix + 1][iy + 1][iz], met0->p[iz]));
106 w[4] =
107 (float) (1e3 * DP2DZ(met0->w[ix][iy][iz + 1], met0->p[iz + 1]));
108 w[5] =
109 (float) (1e3 *
110 DP2DZ(met0->w[ix + 1][iy][iz + 1], met0->p[iz + 1]));
111 w[6] =
112 (float) (1e3 *
113 DP2DZ(met0->w[ix][iy + 1][iz + 1], met0->p[iz + 1]));
114 w[7] =
115 (float) (1e3 *
116 DP2DZ(met0->w[ix + 1][iy + 1][iz + 1], met0->p[iz + 1]));
117
118 /* Collect local wind data... */
119 u[8] = met1->u[ix][iy][iz];
120 u[9] = met1->u[ix + 1][iy][iz];
121 u[10] = met1->u[ix][iy + 1][iz];
122 u[11] = met1->u[ix + 1][iy + 1][iz];
123 u[12] = met1->u[ix][iy][iz + 1];
124 u[13] = met1->u[ix + 1][iy][iz + 1];
125 u[14] = met1->u[ix][iy + 1][iz + 1];
126 u[15] = met1->u[ix + 1][iy + 1][iz + 1];
127
128 v[8] = met1->v[ix][iy][iz];
129 v[9] = met1->v[ix + 1][iy][iz];
130 v[10] = met1->v[ix][iy + 1][iz];
131 v[11] = met1->v[ix + 1][iy + 1][iz];
132 v[12] = met1->v[ix][iy][iz + 1];
133 v[13] = met1->v[ix + 1][iy][iz + 1];
134 v[14] = met1->v[ix][iy + 1][iz + 1];
135 v[15] = met1->v[ix + 1][iy + 1][iz + 1];
136
137 w[8] = (float) (1e3 * DP2DZ(met1->w[ix][iy][iz], met1->p[iz]));
138 w[9] = (float) (1e3 * DP2DZ(met1->w[ix + 1][iy][iz], met1->p[iz]));
139 w[10] = (float) (1e3 * DP2DZ(met1->w[ix][iy + 1][iz], met1->p[iz]));
140 w[11] =
141 (float) (1e3 * DP2DZ(met1->w[ix + 1][iy + 1][iz], met1->p[iz]));
142 w[12] =
143 (float) (1e3 * DP2DZ(met1->w[ix][iy][iz + 1], met1->p[iz + 1]));
144 w[13] =
145 (float) (1e3 *
146 DP2DZ(met1->w[ix + 1][iy][iz + 1], met1->p[iz + 1]));
147 w[14] =
148 (float) (1e3 *
149 DP2DZ(met1->w[ix][iy + 1][iz + 1], met1->p[iz + 1]));
150 w[15] =
151 (float) (1e3 *
152 DP2DZ(met1->w[ix + 1][iy + 1][iz + 1], met1->p[iz + 1]));
153
154 /* Get standard deviations of local wind data... */
155 usig[iz][iy] += stddev(u, 16);
156 vsig[iz][iy] += stddev(v, 16);
157 wsig[iz][iy] += stddev(w, 16);
158 n[iz][iy]++;
159
160 /* Check surface pressure... */
161 if (met0->p[iz] > met0->ps[ix][iy]
162 || met1->p[iz] > met1->ps[ix][iy]) {
163 usig[iz][iy] = NAN;
164 vsig[iz][iy] = NAN;
165 wsig[iz][iy] = NAN;
166 n[iz][iy] = 0;
167 }
168 }
169 }
170
171 /* Create output file... */
172 LOG(1, "Write subgrid data file: %s", argv[2]);
173 if (!(out = fopen(argv[2], "w")))
174 ERRMSG("Cannot create file!");
175
176 /* Write header... */
177 fprintf(out,
178 "# $1 = time [s]\n"
179 "# $2 = altitude [km]\n"
180 "# $3 = longitude [deg]\n"
181 "# $4 = latitude [deg]\n"
182 "# $5 = zonal wind standard deviation [m/s]\n"
183 "# $6 = meridional wind standard deviation [m/s]\n"
184 "# $7 = vertical velocity standard deviation [m/s]\n"
185 "# $8 = number of data points\n");
186
187 /* Write output... */
188 for (int iy = 0; iy < met0->ny - 1; iy++) {
189 fprintf(out, "\n");
190 for (int iz = 0; iz < met0->np - 1; iz++)
191 fprintf(out, "%.2f %g %g %g %g %g %g %d\n",
192 0.5 * (met0->time + met1->time),
193 0.5 * (Z(met0->p[iz]) + Z(met1->p[iz + 1])),
194 0.0, 0.5 * (met0->lat[iy] + met1->lat[iy + 1]),
195 usig[iz][iy] / n[iz][iy], vsig[iz][iy] / n[iz][iy],
196 wsig[iz][iy] / n[iz][iy], n[iz][iy]);
197 }
198
199 /* Close file... */
200 fclose(out);
201
202 /* Free... */
203 free(clim);
204 free(met0);
205 free(met1);
206
207 return EXIT_SUCCESS;
208}
int main(int argc, char *argv[])
Definition: met_subgrid.c:31
float stddev(const float *data, const int n)
Calculates the standard deviation of a set of data.
Definition: mptrac.c:11051
void mptrac_read_clim(const ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
Definition: mptrac.c:5406
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.
Definition: mptrac.c:6357
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:5466
MPTRAC library declarations.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:2043
#define EY
Maximum number of latitudes for meteo data.
Definition: mptrac.h:293
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:1880
#define DP2DZ(dp, p)
Convert a pressure difference to a height difference in the vertical direction.
Definition: mptrac.h:592
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:416
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:1973
#define EP
Maximum number of pressure levels for meteo data.
Definition: mptrac.h:283
Climatological data.
Definition: mptrac.h:3487
Control parameters.
Definition: mptrac.h:2264
Domain decomposition data structure.
Definition: mptrac.h:3720
Meteo data structure.
Definition: mptrac.h:3546
float w[EX][EY][EP]
Vertical velocity [hPa/s].
Definition: mptrac.h:3669
int nx
Number of longitudes.
Definition: mptrac.h:3552
int ny
Number of latitudes.
Definition: mptrac.h:3555
float ps[EX][EY]
Surface pressure [hPa].
Definition: mptrac.h:3585
int np
Number of pressure levels.
Definition: mptrac.h:3558
float u[EX][EY][EP]
Zonal wind [m/s].
Definition: mptrac.h:3663
float v[EX][EY][EP]
Meridional wind [m/s].
Definition: mptrac.h:3666
double time
Time [s].
Definition: mptrac.h:3549
double lat[EY]
Latitudes [deg].
Definition: mptrac.h:3567
double p[EP]
Pressure levels [hPa].
Definition: mptrac.h:3570