MPTRAC
met_zm.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-2024 Forschungszentrum Juelich GmbH
18*/
19
25#include "mptrac.h"
26
27/* ------------------------------------------------------------
28 Dimensions...
29 ------------------------------------------------------------ */
30
32#define NZ 1000
33
35#define NY EY
36
37/* ------------------------------------------------------------
38 Main...
39 ------------------------------------------------------------ */
40
41int main(
42 int argc,
43 char *argv[]) {
44
45 ctl_t ctl;
46
47 clim_t *clim;
48
49 met_t *met;
50
51 FILE *out;
52
53 static double timem[NZ][NY], psm[NZ][NY], tsm[NZ][NY], zsm[NZ][NY],
54 usm[NZ][NY], vsm[NZ][NY], lsmm[NZ][NY], sstm[NZ][NY], pblm[NZ][NY],
55 ptm[NZ][NY], pctm[NZ][NY], pcbm[NZ][NY], clm[NZ][NY], plclm[NZ][NY],
56 plfcm[NZ][NY], pelm[NZ][NY], capem[NZ][NY], cinm[NZ][NY], o3cm[NZ][NY],
57 ttm[NZ][NY], ztm[NZ][NY], tm[NZ][NY], um[NZ][NY], vm[NZ][NY], wm[NZ][NY],
58 h2om[NZ][NY], h2otm[NZ][NY], pvm[NZ][NY], o3m[NZ][NY], lwcm[NZ][NY],
59 rwcm[NZ][NY], iwcm[NZ][NY], swcm[NZ][NY], ccm[NZ][NY], zm[NZ][NY],
60 rhm[NZ][NY], rhicem[NZ][NY], tdewm[NZ][NY], ticem[NZ][NY], tnatm[NZ][NY],
61 hno3m[NZ][NY], ohm[NZ][NY], h2o2m[NZ][NY], ho2m[NZ][NY], o1dm[NZ][NY], z,
62 zt, tt, plev[NZ], ps, ts, zs, us, vs, lsm, sst, pbl, pt, pct, pcb, plcl,
63 plfc, pel, cape, cin, o3c, cl, t, u, v, w, pv, h2o, h2ot, o3, lwc, rwc,
64 iwc, swc, cc, lat, lats[NY], lonm[NZ][NY], cw[3];
65
66 static int np[NZ][NY], npc[NZ][NY], npt[NZ][NY], ny, nz, ci[3];
67
68 /* Allocate... */
69 ALLOC(clim, clim_t, 1);
70 ALLOC(met, met_t, 1);
71
72 /* Check arguments... */
73 if (argc < 4)
74 ERRMSG("Give parameters: <ctl> <zm.tab> <met0> [ <met1> ... ]");
75
76 /* Read control parameters... */
77 read_ctl(argv[1], argc, argv, &ctl);
78 double z0 = scan_ctl(argv[1], argc, argv, "ZM_Z0", -1, "-999", NULL);
79 double z1 = scan_ctl(argv[1], argc, argv, "ZM_Z1", -1, "-999", NULL);
80 double dz = scan_ctl(argv[1], argc, argv, "ZM_DZ", -1, "-999", NULL);
81 double lon0 = scan_ctl(argv[1], argc, argv, "ZM_LON0", -1, "-360", NULL);
82 double lon1 = scan_ctl(argv[1], argc, argv, "ZM_LON1", -1, "360", NULL);
83 double lat0 = scan_ctl(argv[1], argc, argv, "ZM_LAT0", -1, "-90", NULL);
84 double lat1 = scan_ctl(argv[1], argc, argv, "ZM_LAT1", -1, "90", NULL);
85 double dlat = scan_ctl(argv[1], argc, argv, "ZM_DLAT", -1, "-999", NULL);
86
87 /* Read climatological data... */
88 read_clim(&ctl, clim);
89
90 /* Loop over files... */
91 for (int i = 3; i < argc; i++) {
92
93 /* Read meteorological data... */
94 if (!read_met(argv[i], &ctl, clim, met))
95 continue;
96
97 /* Set vertical grid... */
98 if (z0 < 0)
99 z0 = Z(met->p[0]);
100 if (z1 < 0)
101 z1 = Z(met->p[met->np - 1]);
102 nz = 0;
103 if (dz < 0) {
104 for (int iz = 0; iz < met->np; iz++)
105 if (Z(met->p[iz]) >= z0 && Z(met->p[iz]) <= z1) {
106 plev[nz] = met->p[iz];
107 if ((++nz) >= NZ)
108 ERRMSG("Too many pressure levels!");
109 }
110 } else
111 for (z = z0; z <= z1; z += dz) {
112 plev[nz] = P(z);
113 if ((++nz) >= NZ)
114 ERRMSG("Too many pressure levels!");
115 }
116
117 /* Set horizontal grid... */
118 if (dlat <= 0)
119 dlat = fabs(met->lat[1] - met->lat[0]);
120 ny = 0;
121 if (lat0 < -90 && lat1 > 90) {
122 lat0 = gsl_stats_min(met->lat, 1, (size_t) met->ny);
123 lat1 = gsl_stats_max(met->lat, 1, (size_t) met->ny);
124 }
125 for (lat = lat0; lat <= lat1; lat += dlat) {
126 lats[ny] = lat;
127 if ((++ny) >= NY)
128 ERRMSG("Too many latitudes!");
129 }
130
131 /* Average... */
132 for (int ix = 0; ix < met->nx; ix++)
133 if (met->lon[ix] >= lon0 && met->lon[ix] <= lon1)
134 for (int iy = 0; iy < ny; iy++)
135 for (int iz = 0; iz < nz; iz++) {
136
137 /* Interpolate meteo data... */
138 INTPOL_SPACE_ALL(plev[iz], met->lon[ix], lats[iy]);
139
140 /* Averaging... */
141 timem[iz][iy] += met->time;
142 lonm[iz][iy] += met->lon[ix];
143 zm[iz][iy] += z;
144 tm[iz][iy] += t;
145 um[iz][iy] += u;
146 vm[iz][iy] += v;
147 wm[iz][iy] += w;
148 pvm[iz][iy] += pv;
149 h2om[iz][iy] += h2o;
150 o3m[iz][iy] += o3;
151 lwcm[iz][iy] += lwc;
152 rwcm[iz][iy] += rwc;
153 iwcm[iz][iy] += iwc;
154 swcm[iz][iy] += swc;
155 ccm[iz][iy] += cc;
156 psm[iz][iy] += ps;
157 tsm[iz][iy] += ts;
158 zsm[iz][iy] += zs;
159 usm[iz][iy] += us;
160 vsm[iz][iy] += vs;
161 lsmm[iz][iy] += lsm;
162 sstm[iz][iy] += sst;
163 pblm[iz][iy] += pbl;
164 pctm[iz][iy] += pct;
165 pcbm[iz][iy] += pcb;
166 clm[iz][iy] += cl;
167 if (isfinite(plfc) && isfinite(pel) && cape >= ctl.conv_cape
168 && (ctl.conv_cin <= 0 || cin < ctl.conv_cin)) {
169 plclm[iz][iy] += plcl;
170 plfcm[iz][iy] += plfc;
171 pelm[iz][iy] += pel;
172 capem[iz][iy] += cape;
173 cinm[iz][iy] += cin;
174 npc[iz][iy]++;
175 }
176 if (isfinite(pt)) {
177 ptm[iz][iy] += pt;
178 ztm[iz][iy] += zt;
179 ttm[iz][iy] += tt;
180 h2otm[iz][iy] += h2ot;
181 npt[iz][iy]++;
182 }
183 o3cm[iz][iy] += o3c;
184 rhm[iz][iy] += RH(plev[iz], t, h2o);
185 rhicem[iz][iy] += RHICE(plev[iz], t, h2o);
186 tdewm[iz][iy] += TDEW(plev[iz], h2o);
187 ticem[iz][iy] += TICE(plev[iz], h2o);
188 hno3m[iz][iy] +=
189 clim_zm(&clim->hno3, met->time, lats[iy], plev[iz]);
190 tnatm[iz][iy] +=
191 nat_temperature(plev[iz], h2o,
192 clim_zm(&clim->hno3, met->time, lats[iy],
193 plev[iz]));
194 ohm[iz][iy] +=
195 clim_oh(&ctl, clim, met->time, met->lon[ix], lats[iy],
196 plev[iz]);
197 h2o2m[iz][iy]
198 += clim_zm(&clim->h2o2, met->time, lats[iy], plev[iz]);
199 ho2m[iz][iy]
200 += clim_zm(&clim->ho2, met->time, lats[iy], plev[iz]);
201 o1dm[iz][iy]
202 += clim_zm(&clim->o1d, met->time, lats[iy], plev[iz]);
203 np[iz][iy]++;
204 }
205 }
206
207 /* Create output file... */
208 LOG(1, "Write meteorological data file: %s", argv[2]);
209 if (!(out = fopen(argv[2], "w")))
210 ERRMSG("Cannot create file!");
211
212 /* Write header... */
213 MET_HEADER;
214
215 /* Write data... */
216 for (int iz = 0; iz < nz; iz++) {
217 fprintf(out, "\n");
218 for (int iy = 0; iy < ny; iy++)
219 fprintf(out,
220 "%.2f %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
221 " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g"
222 " %g %g %g %g %g %g %g %g %g %g %d %d %d\n",
223 timem[iz][iy] / np[iz][iy], Z(plev[iz]),
224 lonm[iz][iy] / np[iz][iy], lats[iy],
225 plev[iz], tm[iz][iy] / np[iz][iy], um[iz][iy] / np[iz][iy],
226 vm[iz][iy] / np[iz][iy], wm[iz][iy] / np[iz][iy],
227 h2om[iz][iy] / np[iz][iy], o3m[iz][iy] / np[iz][iy],
228 zm[iz][iy] / np[iz][iy], pvm[iz][iy] / np[iz][iy],
229 psm[iz][iy] / np[iz][iy], tsm[iz][iy] / np[iz][iy],
230 zsm[iz][iy] / np[iz][iy], usm[iz][iy] / np[iz][iy],
231 vsm[iz][iy] / np[iz][iy], lsmm[iz][iy] / np[iz][iy],
232 sstm[iz][iy] / np[iz][iy], ptm[iz][iy] / npt[iz][iy],
233 ztm[iz][iy] / npt[iz][iy], ttm[iz][iy] / npt[iz][iy],
234 h2otm[iz][iy] / npt[iz][iy], lwcm[iz][iy] / np[iz][iy],
235 rwcm[iz][iy] / np[iz][iy], iwcm[iz][iy] / np[iz][iy],
236 swcm[iz][iy] / np[iz][iy], ccm[iz][iy] / np[iz][iy],
237 clm[iz][iy] / np[iz][iy], pctm[iz][iy] / np[iz][iy],
238 pcbm[iz][iy] / np[iz][iy], plclm[iz][iy] / npc[iz][iy],
239 plfcm[iz][iy] / npc[iz][iy], pelm[iz][iy] / npc[iz][iy],
240 capem[iz][iy] / npc[iz][iy], cinm[iz][iy] / npc[iz][iy],
241 rhm[iz][iy] / np[iz][iy], rhicem[iz][iy] / np[iz][iy],
242 tdewm[iz][iy] / np[iz][iy], ticem[iz][iy] / np[iz][iy],
243 tnatm[iz][iy] / np[iz][iy], hno3m[iz][iy] / np[iz][iy],
244 ohm[iz][iy] / np[iz][iy], h2o2m[iz][iy] / np[iz][iy],
245 ho2m[iz][iy] / np[iz][iy], o1dm[iz][iy] / np[iz][iy],
246 pblm[iz][iy] / np[iz][iy], o3cm[iz][iy] / np[iz][iy],
247 np[iz][iy], npt[iz][iy], npc[iz][iy]);
248 }
249
250 /* Close file... */
251 fclose(out);
252
253 /* Free... */
254 free(clim);
255 free(met);
256
257 return EXIT_SUCCESS;
258}
main
int main(int argc, char *argv[])
Definition: met_zm.c:41
NY
#define NY
Maximum number of latitudes.
Definition: met_zm.c:35
NZ
#define NZ
Maximum number of altitudes.
Definition: met_zm.c:32
clim_zm
double clim_zm(const clim_zm_t *zm, const double t, const double lat, const double p)
Interpolates monthly mean zonal mean climatological variables.
Definition: mptrac.c:401
nat_temperature
double nat_temperature(const double p, const double h2o, const double hno3)
Calculates the nitric acid trihydrate (NAT) temperature.
Definition: mptrac.c:4191
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:8098
clim_oh
double clim_oh(const ctl_t *ctl, const clim_t *clim, const double t, const double lon, const double lat, const double p)
Calculates the hydroxyl radical (OH) concentration from climatology data, with an optional diurnal co...
Definition: mptrac.c:89
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:4805
read_met
int read_met(const char *filename, ctl_t *ctl, clim_t *clim, met_t *met)
Reads meteorological data from a file, supporting multiple formats and MPI broadcasting.
Definition: mptrac.c:5655
read_clim
void read_clim(const ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
Definition: mptrac.c:4473
mptrac.h
MPTRAC library declarations.
INTPOL_SPACE_ALL
#define INTPOL_SPACE_ALL(p, lon, lat)
Interpolate multiple meteorological variables in space.
Definition: mptrac.h:726
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:1901
Z
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:1726
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1289
MET_HEADER
#define MET_HEADER
Write header for meteorological data file.
Definition: mptrac.h:888
TICE
#define TICE(p, h2o)
Calculate frost point temperature (WMO, 2018).
Definition: mptrac.h:1605
RHICE
#define RHICE(p, t, h2o)
Compute relative humidity over ice.
Definition: mptrac.h:1441
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:347
RH
#define RH(p, t, h2o)
Compute relative humidity over water.
Definition: mptrac.h:1411
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:1831
TDEW
#define TDEW(p, h2o)
Calculate dew point temperature.
Definition: mptrac.h:1580
clim_t
Climatological data.
Definition: mptrac.h:3282
clim_t::ho2
clim_zm_t ho2
HO2 zonal means.
Definition: mptrac.h:3312
clim_t::hno3
clim_zm_t hno3
HNO3 zonal means.
Definition: mptrac.h:3303
clim_t::o1d
clim_zm_t o1d
O(1D) zonal means.
Definition: mptrac.h:3315
clim_t::h2o2
clim_zm_t h2o2
H2O2 zonal means.
Definition: mptrac.h:3309
ctl_t
Control parameters.
Definition: mptrac.h:2155
ctl_t::conv_cape
double conv_cape
CAPE threshold for convection module [J/kg].
Definition: mptrac.h:2663
ctl_t::conv_cin
double conv_cin
CIN threshold for convection module [J/kg].
Definition: mptrac.h:2666
met_t
Meteo data structure.
Definition: mptrac.h:3341
met_t::nx
int nx
Number of longitudes.
Definition: mptrac.h:3347
met_t::ny
int ny
Number of latitudes.
Definition: mptrac.h:3350
met_t::np
int np
Number of pressure levels.
Definition: mptrac.h:3353
met_t::lon
double lon[EX]
Longitude [deg].
Definition: mptrac.h:3359
met_t::time
double time
Time [s].
Definition: mptrac.h:3344
met_t::lat
double lat[EY]
Latitude [deg].
Definition: mptrac.h:3362
met_t::p
double p[EP]
Pressure levels [hPa].
Definition: mptrac.h:3365
Generated on Thu Nov 21 2024 13:53:30 for MPTRAC by doxygen 1.9.4