41 {
42
44
46
48
50
53
54 static int init, nx, ny, nt, ncid, varid, dims[3];
55
56 static size_t count[10], start[10];
57
58
62
63
65
66
67 if (argc < 4)
68 ERRMSG(
"Missing or invalid command-line arguments.\n\n"
69 "Usage: tropo <ctl> <tropo.nc> <met0> [<met1> ...]\n\n"
70 "Use -h for full help.");
71
72
74 double lon0 =
scan_ctl(argv[1], argc, argv,
"TROPO_LON0", -1,
"-180", NULL);
75 double lon1 =
scan_ctl(argv[1], argc, argv,
"TROPO_LON1", -1,
"180", NULL);
76 double dlon =
scan_ctl(argv[1], argc, argv,
"TROPO_DLON", -1,
"-999", NULL);
77 double lat0 =
scan_ctl(argv[1], argc, argv,
"TROPO_LAT0", -1,
"-90", NULL);
78 double lat1 =
scan_ctl(argv[1], argc, argv,
"TROPO_LAT1", -1,
"90", NULL);
79 double dlat =
scan_ctl(argv[1], argc, argv,
"TROPO_DLAT", -1,
"-999", NULL);
80 int h2o = (int)
scan_ctl(argv[1], argc, argv,
"TROPO_H2O", -1,
"1", NULL);
81 int o3 = (int)
scan_ctl(argv[1], argc, argv,
"TROPO_O3", -1,
"1", NULL);
82
83
85
86
87 for (int i = 3; i < argc; i++) {
88
89
91
92
94 continue;
95
96
97 if (!init) {
98 init = 1;
99
100
101 if (dlon <= 0)
102 dlon = fabs(met->
lon[1] - met->
lon[0]);
103 if (dlat <= 0)
104 dlat = fabs(met->
lat[1] - met->
lat[0]);
105 if (lon0 < -360 && lon1 > 360) {
106 lon0 = gsl_stats_min(met->
lon, 1, (
size_t) met->
nx);
107 lon1 = gsl_stats_max(met->
lon, 1, (
size_t) met->
nx);
108 }
109 nx = ny = 0;
110 for (lon = lon0; lon <= lon1 + 0.001; lon += dlon) {
111 lons[nx] = round(lon * 1e3) / 1e3;
113 ERRMSG(
"Too many longitudes!");
114 }
115 if (lat0 < -90 && lat1 > 90) {
116 lat0 = gsl_stats_min(met->
lat, 1, (
size_t) met->
ny);
117 lat1 = gsl_stats_max(met->
lat, 1, (
size_t) met->
ny);
118 }
119 for (lat = lat0; lat <= lat1 + 0.001; lat += dlat) {
120 lats[ny] = round(lat * 1e3) / 1e3;
122 ERRMSG(
"Too many latitudes!");
123 }
124
125
126 LOG(1,
"Write tropopause data file: %s", argv[2]);
127 NC(nc_create(argv[2], NC_NETCDF4, &ncid));
128
129
130 NC(nc_def_dim(ncid,
"time", (
size_t) NC_UNLIMITED, &dims[0]));
131 NC(nc_def_dim(ncid,
"lat", (
size_t) ny, &dims[1]));
132 NC(nc_def_dim(ncid,
"lon", (
size_t) nx, &dims[2]));
133
134
135 NC_DEF_VAR(
"time", NC_DOUBLE, 1, &dims[0],
"time",
136 "seconds since 2000-01-01 00:00:00 UTC", 0, 0);
137 NC_DEF_VAR(
"lat", NC_DOUBLE, 1, &dims[1],
"latitude",
"degrees_north",
138 0, 0);
139 NC_DEF_VAR(
"lon", NC_DOUBLE, 1, &dims[2],
"longitude",
"degrees_east",
140 0, 0);
141
142 NC_DEF_VAR(
"clp_z", NC_FLOAT, 3, &dims[0],
"cold point height",
"km", 0,
143 0);
144 NC_DEF_VAR(
"clp_p", NC_FLOAT, 3, &dims[0],
"cold point pressure",
"hPa",
145 0, 0);
146 NC_DEF_VAR(
"clp_t", NC_FLOAT, 3, &dims[0],
"cold point temperature",
147 "K", 0, 0);
148 if (h2o)
149 NC_DEF_VAR(
"clp_q", NC_FLOAT, 3, &dims[0],
"cold point water vapor",
150 "ppv", 0, 0);
151 if (o3)
152 NC_DEF_VAR(
"clp_o3", NC_FLOAT, 3, &dims[0],
"cold point ozone",
153 "ppv", 0, 0);
154
156 "dynamical tropopause height", "km", 0, 0);
158 "dynamical tropopause pressure", "hPa", 0, 0);
160 "dynamical tropopause temperature", "K", 0, 0);
161 if (h2o)
163 "dynamical tropopause water vapor", "ppv", 0, 0);
164 if (o3)
166 "dynamical tropopause ozone", "ppv", 0, 0);
167
168 NC_DEF_VAR(
"wmo_1st_z", NC_FLOAT, 3, &dims[0],
169 "WMO 1st tropopause height", "km", 0, 0);
170 NC_DEF_VAR(
"wmo_1st_p", NC_FLOAT, 3, &dims[0],
171 "WMO 1st tropopause pressure", "hPa", 0, 0);
172 NC_DEF_VAR(
"wmo_1st_t", NC_FLOAT, 3, &dims[0],
173 "WMO 1st tropopause temperature", "K", 0, 0);
174 if (h2o)
175 NC_DEF_VAR(
"wmo_1st_q", NC_FLOAT, 3, &dims[0],
176 "WMO 1st tropopause water vapor", "ppv", 0, 0);
177 if (o3)
178 NC_DEF_VAR(
"wmo_1st_o3", NC_FLOAT, 3, &dims[0],
179 "WMO 1st tropopause ozone", "ppv", 0, 0);
180
181 NC_DEF_VAR(
"wmo_2nd_z", NC_FLOAT, 3, &dims[0],
182 "WMO 2nd tropopause height", "km", 0, 0);
183 NC_DEF_VAR(
"wmo_2nd_p", NC_FLOAT, 3, &dims[0],
184 "WMO 2nd tropopause pressure", "hPa", 0, 0);
185 NC_DEF_VAR(
"wmo_2nd_t", NC_FLOAT, 3, &dims[0],
186 "WMO 2nd tropopause temperature", "K", 0, 0);
187 if (h2o)
188 NC_DEF_VAR(
"wmo_2nd_q", NC_FLOAT, 3, &dims[0],
189 "WMO 2nd tropopause water vapor", "ppv", 0, 0);
190 if (o3)
191 NC_DEF_VAR(
"wmo_2nd_o3", NC_FLOAT, 3, &dims[0],
192 "WMO 2nd tropopause ozone", "ppv", 0, 0);
193
194 NC_DEF_VAR(
"ps", NC_FLOAT, 3, &dims[0],
"surface pressure",
"hPa", 0,
195 0);
196 NC_DEF_VAR(
"zs", NC_FLOAT, 3, &dims[0],
"surface height",
"km", 0, 0);
197
198
200
201
204 }
205
206
207 start[0] = (size_t) nt;
208 count[0] = 1;
209 start[1] = 0;
210 count[1] = (size_t) ny;
211 start[2] = 0;
212 count[2] = (size_t) nx;
214
215
216 get_tropo(2, &ctl, clim, met, lons, nx, lats, ny, pt, zt, tt, qt, o3t, ps,
217 zs);
221 if (h2o)
223 if (o3)
225
226
227 get_tropo(5, &ctl, clim, met, lons, nx, lats, ny, pt, zt, tt, qt, o3t, ps,
228 zs);
232 if (h2o)
234 if (o3)
236
237
238 get_tropo(3, &ctl, clim, met, lons, nx, lats, ny, pt, zt, tt, qt, o3t, ps,
239 zs);
243 if (h2o)
245 if (o3)
247
248
249 get_tropo(4, &ctl, clim, met, lons, nx, lats, ny, pt, zt, tt, qt, o3t, ps,
250 zs);
254 if (h2o)
256 if (o3)
258
259
262
263
264 nt++;
265 }
266
267
269
270
271 free(clim);
272 free(met);
273 free(dd);
274
275 return EXIT_SUCCESS;
276}
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 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.
void get_tropo(const int met_tropo, ctl_t *ctl, const clim_t *clim, met_t *met, const double *lons, const int nx, const double *lats, const int ny, double *pt, double *zt, double *tt, double *qt, double *o3t, double *ps, double *zs)
Calculate tropopause data.
#define NC(cmd)
Execute a NetCDF command and check for errors.
#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 EX
Maximum number of longitudes for meteo data.
#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 NC_DEF_VAR(varname, type, ndims, dims, long_name, units, level, quant)
Define a NetCDF variable with attributes.
#define NC_PUT_DOUBLE(varname, ptr, hyperslab)
Write double precision data to a NetCDF variable.
int met_tropo
Tropopause definition (0=none, 1=clim, 2=cold point, 3=WMO_1st, 4=WMO_2nd, 5=dynamical).
Domain decomposition data structure.
int nx
Number of longitudes.
int ny
Number of latitudes.
double lon[EX]
Longitudes [deg].
double lat[EY]
Latitudes [deg].