MPTRAC
Macros | Functions
tropo_zm.c File Reference

Extract zonal mean of tropopause data set. More...

#include "mptrac.h"

Go to the source code of this file.

Macros

#define NT   744
 Maximum number of time steps. More...
 

Functions

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

Detailed Description

Extract zonal mean of tropopause data set.

Definition in file tropo_zm.c.

Macro Definition Documentation

◆ NT

#define NT   744

Maximum number of time steps.

Definition at line 32 of file tropo_zm.c.

Function Documentation

◆ main()

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

Definition at line 38 of file tropo_zm.c.

40 {
41
42 ctl_t ctl;
43
44 static FILE *out;
45
46 static char varname[LEN];
47
48 static double time0, lons[EX], lats[EY], zm[EY], zs[EY], pm[EY],
49 ps[EY], tm[EY], ts[EY], qm[EY], qs[EY], o3m[EY], o3s[EY];
50
51 static float help[EX * EY], tropo_z0[EX][EY], tropo_p0[EX][EY],
52 tropo_t0[EX][EY], tropo_q0[EX][EY], tropo_o30[EX][EY];
53
54 static int ncid, varid, varid_z, varid_p, varid_t, varid_q, varid_o3, h2o,
55 o3, n[EY], nt[EY], init, ntime, nlon, nlat, ilon, ilat;
56
57 static size_t count[10], start[10];
58
59 /* Check arguments... */
60 if (argc < 5)
61 ERRMSG("Give parameters: <ctl> <zm.tab> <var> <tropo.nc>");
62
63 /* Read control parameters... */
64 read_ctl(argv[1], argc, argv, &ctl);
65
66 /* Loop over tropopause files... */
67 for (int iarg = 4; iarg < argc; iarg++) {
68
69 /* Open tropopause file... */
70 LOG(1, "Read tropopause data: %s", argv[iarg]);
71 if (nc_open(argv[iarg], NC_NOWRITE, &ncid) != NC_NOERR)
72 ERRMSG("Cannot open file!");
73
74 /* Get dimensions... */
75 NC_INQ_DIM("time", &ntime, 1, NT);
76 NC_INQ_DIM("lat", &nlat, 1, EY);
77 NC_INQ_DIM("lon", &nlon, 1, EX);
78
79 /* Read coordinates... */
80 NC_GET_DOUBLE("lat", lats, 1);
81 NC_GET_DOUBLE("lon", lons, 1);
82
83 /* Get variable indices... */
84 sprintf(varname, "%s_z", argv[3]);
85 NC(nc_inq_varid(ncid, varname, &varid_z));
86 sprintf(varname, "%s_p", argv[3]);
87 NC(nc_inq_varid(ncid, varname, &varid_p));
88 sprintf(varname, "%s_t", argv[3]);
89 NC(nc_inq_varid(ncid, varname, &varid_t));
90 sprintf(varname, "%s_q", argv[3]);
91 h2o = (nc_inq_varid(ncid, varname, &varid_q) == NC_NOERR);
92 sprintf(varname, "%s_o3", argv[3]);
93 o3 = (nc_inq_varid(ncid, varname, &varid_o3) == NC_NOERR);
94
95 /* Set dimensions... */
96 count[0] = 1;
97 count[1] = (size_t) nlat;
98 count[2] = (size_t) nlon;
99
100 /* Loop over time steps... */
101 for (int it = 0; it < ntime; it++) {
102
103 /* Get time from filename... */
104 if (!init) {
105 init = 1;
106 time0 = time_from_filename(argv[iarg], 13);
107 }
108
109 /* Read data... */
110 start[0] = (size_t) it;
111 NC(nc_get_vara_float(ncid, varid_z, start, count, help));
112 for (ilon = 0; ilon < nlon; ilon++)
113 for (ilat = 0; ilat < nlat; ilat++)
114 tropo_z0[ilon][ilat] = help[ilat * nlon + ilon];
115 NC(nc_get_vara_float(ncid, varid_p, start, count, help));
116 for (ilon = 0; ilon < nlon; ilon++)
117 for (ilat = 0; ilat < nlat; ilat++)
118 tropo_p0[ilon][ilat] = help[ilat * nlon + ilon];
119 NC(nc_get_vara_float(ncid, varid_t, start, count, help));
120 for (ilon = 0; ilon < nlon; ilon++)
121 for (ilat = 0; ilat < nlat; ilat++)
122 tropo_t0[ilon][ilat] = help[ilat * nlon + ilon];
123 if (h2o) {
124 NC(nc_get_vara_float(ncid, varid_q, start, count, help));
125 for (ilon = 0; ilon < nlon; ilon++)
126 for (ilat = 0; ilat < nlat; ilat++)
127 tropo_q0[ilon][ilat] = help[ilat * nlon + ilon];
128 } else
129 for (ilon = 0; ilon < nlon; ilon++)
130 for (ilat = 0; ilat < nlat; ilat++)
131 tropo_q0[ilon][ilat] = NAN;
132 if (o3) {
133 NC(nc_get_vara_float(ncid, varid_o3, start, count, help));
134 for (ilon = 0; ilon < nlon; ilon++)
135 for (ilat = 0; ilat < nlat; ilat++)
136 tropo_o30[ilon][ilat] = help[ilat * nlon + ilon];
137 } else
138 for (ilon = 0; ilon < nlon; ilon++)
139 for (ilat = 0; ilat < nlat; ilat++)
140 tropo_o30[ilon][ilat] = NAN;
141
142 /* Averaging... */
143 for (ilat = 0; ilat < nlat; ilat++)
144 for (ilon = 0; ilon < nlon; ilon++) {
145 nt[ilat]++;
146 if (isfinite(tropo_z0[ilon][ilat])
147 && isfinite(tropo_p0[ilon][ilat])
148 && isfinite(tropo_t0[ilon][ilat])
149 && (!h2o || isfinite(tropo_q0[ilon][ilat]))
150 && (!o3 || isfinite(tropo_o30[ilon][ilat]))) {
151 zm[ilat] += tropo_z0[ilon][ilat];
152 zs[ilat] += SQR(tropo_z0[ilon][ilat]);
153 pm[ilat] += tropo_p0[ilon][ilat];
154 ps[ilat] += SQR(tropo_p0[ilon][ilat]);
155 tm[ilat] += tropo_t0[ilon][ilat];
156 ts[ilat] += SQR(tropo_t0[ilon][ilat]);
157 qm[ilat] += tropo_q0[ilon][ilat];
158 qs[ilat] += SQR(tropo_q0[ilon][ilat]);
159 o3m[ilat] += tropo_o30[ilon][ilat];
160 o3s[ilat] += SQR(tropo_o30[ilon][ilat]);
161 n[ilat]++;
162 }
163 }
164 }
165
166 /* Close files... */
167 NC(nc_close(ncid));
168 }
169
170 /* Normalize... */
171 for (ilat = 0; ilat < nlat; ilat++)
172 if (n[ilat] > 0) {
173 zm[ilat] /= n[ilat];
174 pm[ilat] /= n[ilat];
175 tm[ilat] /= n[ilat];
176 qm[ilat] /= n[ilat];
177 o3m[ilat] /= n[ilat];
178 double aux = zs[ilat] / n[ilat] - SQR(zm[ilat]);
179 zs[ilat] = aux > 0 ? sqrt(aux) : 0.0;
180 aux = ps[ilat] / n[ilat] - SQR(pm[ilat]);
181 ps[ilat] = aux > 0 ? sqrt(aux) : 0.0;
182 aux = ts[ilat] / n[ilat] - SQR(tm[ilat]);
183 ts[ilat] = aux > 0 ? sqrt(aux) : 0.0;
184 aux = qs[ilat] / n[ilat] - SQR(qm[ilat]);
185 qs[ilat] = aux > 0 ? sqrt(aux) : 0.0;
186 aux = o3s[ilat] / n[ilat] - SQR(o3m[ilat]);
187 o3s[ilat] = aux > 0 ? sqrt(aux) : 0.0;
188 }
189
190 /* Create file... */
191 LOG(1, "Write tropopause zonal mean data: %s", argv[2]);
192 if (!(out = fopen(argv[2], "w")))
193 ERRMSG("Cannot create file!");
194
195 /* Write header... */
196 fprintf(out,
197 "# $1 = time [s]\n"
198 "# $2 = latitude [deg]\n"
199 "# $3 = tropopause height (mean) [km]\n"
200 "# $4 = tropopause pressure (mean) [hPa]\n"
201 "# $5 = tropopause temperature (mean) [K]\n"
202 "# $6 = tropopause water vapor (mean) [ppv]\n"
203 "# $7 = tropopause ozone (mean) [ppv]\n"
204 "# $8 = tropopause height (sigma) [km]\n"
205 "# $9 = tropopause pressure (sigma) [hPa]\n"
206 "# $10 = tropopause temperature (sigma) [K]\n"
207 "# $11 = tropopause water vapor (sigma) [ppv]\n"
208 "# $12 = tropopause ozone (sigma) [ppv]\n"
209 "# $13 = number of data points\n"
210 "# $14 = occurrence frequency [%%]\n\n");
211
212 /* Write output... */
213 for (ilat = 0; ilat < nlat; ilat++)
214 fprintf(out, "%.2f %g %g %g %g %g %g %g %g %g %g %g %d %g\n",
215 time0, lats[ilat], zm[ilat], pm[ilat], tm[ilat], qm[ilat],
216 o3m[ilat], zs[ilat], ps[ilat], ts[ilat], qs[ilat], o3s[ilat],
217 n[ilat], 100. * n[ilat] / nt[ilat]);
218
219 /* Close files... */
220 fclose(out);
221
222 return EXIT_SUCCESS;
223}
time_from_filename
double time_from_filename(const char *filename, int offset)
Extracts and converts a timestamp from a filename to Julian seconds.
Definition: mptrac.c:8258
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
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: mptrac.h:236
NC
#define NC(cmd)
Execute a NetCDF command and check for errors.
Definition: mptrac.h:981
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:1881
EY
#define EY
Maximum number of latitudes for meteo data.
Definition: mptrac.h:266
EX
#define EX
Maximum number of longitudes for meteo data.
Definition: mptrac.h:261
SQR
#define SQR(x)
Compute the square of a value.
Definition: mptrac.h:1522
NC_INQ_DIM
#define NC_INQ_DIM(dimname, ptr, min, max)
Inquire the length of a dimension in a NetCDF file.
Definition: mptrac.h:1054
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:1811
NC_GET_DOUBLE
#define NC_GET_DOUBLE(varname, ptr, force)
Retrieve a double-precision variable from a NetCDF file.
Definition: mptrac.h:1026
ctl_t
Control parameters.
Definition: mptrac.h:2135
NT
#define NT
Maximum number of time steps.
Definition: tropo_zm.c:32
Here is the call graph for this function:
Generated on Mon Sep 16 2024 17:36:26 for MPTRAC by doxygen 1.9.4