46 {
47
49
51
53
54 FILE *out;
55
56 static double p2[1000], t[1000], t2[1000], z[1000], z2[1000], lat_mean,
57 z_mean;
58
59 static int hist_max[1000], hist_min[1000], hist_mean[1000], hist_sig[1000],
60 nhist_max, nhist_min, nhist_mean, nhist_sig, np;
61
62
65
66
67 if (argc < 4)
68 ERRMSG(
"Give parameters: <ctl> <lapse.tab> <met0> [ <met1> ... ]");
69
70
72 const int dz =
73 (int)
scan_ctl(argv[1], argc, argv,
"LAPSE_DZ", -1,
"20", NULL);
74 const double lat0 =
75 (int)
scan_ctl(argv[1], argc, argv,
"LAPSE_LAT0", -1,
"-90", NULL);
76 const double lat1 =
77 (int)
scan_ctl(argv[1], argc, argv,
"LAPSE_LAT1", -1,
"90", NULL);
78 const double z0 =
79 (int)
scan_ctl(argv[1], argc, argv,
"LAPSE_Z0", -1,
"0", NULL);
80 const double z1 =
81 (int)
scan_ctl(argv[1], argc, argv,
"LAPSE_Z1", -1,
"100", NULL);
82 const int intpol =
83 (int)
scan_ctl(argv[1], argc, argv,
"LAPSE_INTPOL", -1,
"1", NULL);
84
85
87
88
89 for (int i = 3; i < argc; i++) {
90
91
92 if (!
read_met(argv[i], &ctl, clim, met))
93 continue;
94
95
96 for (
int iz = 0; iz < met->
np; iz++)
97 z[iz] =
Z(met->
p[iz]);
98 for (int iz = 0; iz <= 250; iz++) {
99 z2[iz] = 0.0 + 0.1 * iz;
101 }
102
103
104 for (
int ix = 0; ix < met->
nx; ix++)
105 for (
int iy = 0; iy < met->
ny; iy++) {
106
107
108 if (met->
lat[iy] < lat0 || met->
lat[iy] > lat1)
109 continue;
110
111
112 for (
int iz = 0; iz < met->
np; iz++)
113 t[iz] = met->
t[ix][iy][iz];
114 if (intpol == 1)
116 else
117 for (int iz = 0; iz <= 250; iz++) {
119 t2[iz] =
LIN(z[idx], t[idx], z[idx + 1], t[idx + 1], z2[iz]);
120 }
121
122
123 for (int iz = 0; iz <= 250; iz++) {
124
125
126 if (z2[iz] < z0 || z2[iz] > z1)
127 continue;
128
129
130 if (p2[iz] > met->
ps[ix][iy])
131 continue;
132
133
134 lat_mean += met->
lat[iy];
135 z_mean += z2[iz];
136 np++;
137
138
139 int nlapse = 0;
140 double lapse_max = -1e99, lapse_min = 1e99, lapse_mean =
141 0, lapse_sig = 0;
142 for (int iz2 = iz + 1; iz2 <= iz + dz; iz2++) {
143 lapse_max =
144 MAX(
LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]), lapse_max);
145 lapse_min =
146 MIN(
LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]), lapse_min);
147 lapse_mean +=
LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]);
148 lapse_sig +=
SQR(
LAPSE(p2[iz], t2[iz], p2[iz2], t2[iz2]));
149 nlapse++;
150 }
151 lapse_mean /= nlapse;
152 lapse_sig = sqrt(
MAX(lapse_sig / nlapse -
SQR(lapse_mean), 0));
153
154
156 if (idx >= 0 && idx <
IDXMAX) {
157 hist_max[idx]++;
158 nhist_max++;
159 }
160
162 if (idx >= 0 && idx <
IDXMAX) {
163 hist_min[idx]++;
164 nhist_min++;
165 }
166
168 if (idx >= 0 && idx <
IDXMAX) {
169 hist_mean[idx]++;
170 nhist_mean++;
171 }
172
174 if (idx >= 0 && idx <
IDXMAX) {
175 hist_sig[idx]++;
176 nhist_sig++;
177 }
178 }
179 }
180 }
181
182
183 LOG(1,
"Write lapse rate data: %s", argv[2]);
184 if (!(out = fopen(argv[2], "w")))
185 ERRMSG(
"Cannot create file!");
186
187
188 fprintf(out,
189 "# $1 = mean altitude [km]\n"
190 "# $2 = mean latitude [deg]\n"
191 "# $3 = lapse rate [K/km]\n"
192 "# $4 = counts of maxima per bin\n"
193 "# $5 = total number of maxima\n"
194 "# $6 = normalized frequency of maxima\n"
195 "# $7 = counts of minima per bin\n"
196 "# $8 = total number of minima\n"
197 "# $9 = normalized frequency of minima\n"
198 "# $10 = counts of means per bin\n"
199 "# $11 = total number of means\n"
200 "# $12 = normalized frequency of means\n"
201 "# $13 = counts of sigmas per bin\n"
202 "# $14 = total number of sigmas\n"
203 "# $15 = normalized frequency of sigmas\n\n");
204
205
206 double nmax_max = 0, nmax_min = 0, nmax_mean = 0, nmax_sig = 0;
207 for (
int idx = 0; idx <
IDXMAX; idx++) {
208 nmax_max =
MAX(hist_max[idx], nmax_max);
209 nmax_min =
MAX(hist_min[idx], nmax_min);
210 nmax_mean =
MAX(hist_mean[idx], nmax_mean);
211 nmax_sig =
MAX(hist_sig[idx], nmax_sig);
212 }
213 for (
int idx = 0; idx <
IDXMAX; idx++)
214 fprintf(out,
215 "%g %g %g %d %d %g %d %d %g %d %d %g %d %d %g\n",
217 hist_max[idx], nhist_max,
218 (double) hist_max[idx] / (double) nmax_max, hist_min[idx],
219 nhist_min, (double) hist_min[idx] / (double) nmax_min,
220 hist_mean[idx], nhist_mean,
221 (double) hist_mean[idx] / (double) nmax_mean, hist_sig[idx],
222 nhist_sig, (double) hist_sig[idx] / (double) nmax_sig);
223
224
225 fclose(out);
226
227
228 free(clim);
229 free(met);
230
231 return EXIT_SUCCESS;
232}
#define LAPSEMIN
Lapse rate minimum [K/km.
#define IDXMAX
Maximum number of histogram bins.
#define DLAPSE
Lapse rate bin size [K/km].
int read_met(const char *filename, const ctl_t *ctl, const clim_t *clim, met_t *met)
Reads meteorological data from a file, supporting multiple formats and MPI broadcasting.
int locate_irr(const double *xx, const int n, const double x)
Locate the index of the interval containing a given value in a sorted array.
void spline(const double *x, const double *y, const int n, const double *x2, double *y2, const int n2, const int method)
Performs spline interpolation or linear interpolation.
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 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 read_clim(const ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
#define LAPSE(p1, t1, p2, t2)
Calculate lapse rate.
#define MIN(a, b)
Macro to determine the minimum of two values.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
#define Z(p)
Convert pressure to altitude.
#define P(z)
Compute pressure at given altitude.
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
#define SQR(x)
Compute the square of a value.
#define LOG(level,...)
Print a log message with a specified logging level.
#define LIN(x0, y0, x1, y1, x)
Linear interpolation.
#define MAX(a, b)
Macro to determine the maximum of two values.
int met_tropo_spline
Tropopause interpolation method (0=linear, 1=spline).
int nx
Number of longitudes.
int ny
Number of latitudes.
float ps[EX][EY]
Surface pressure [hPa].
int np
Number of pressure levels.
float t[EX][EY][EP]
Temperature [K].
double lat[EY]
Latitude [deg].
double p[EP]
Pressure levels [hPa].