AIRS Code Collection
Macros | Functions
hurricane.c File Reference

Analyse gravity wave data for tropical cyclones. More...

#include "libairs.h"

Go to the source code of this file.

Macros

#define NSTORM   9000
 
#define NTIME   140
 

Functions

int get_storm_pos (int nobs, double time_wmo[NTIME], double lon_wmo[NTIME], double lat_wmo[NTIME], double wind_wmo[NTIME], double pres_wmo[NTIME], double t, int dt, int st, double x[3], double *wind, double *dwind, double *pres, double *dpres)
 
void read_var (int ncid, const char varname[], size_t nstorm, int nobs[NSTORM], double x[NSTORM][NTIME])
 
int main (int argc, char *argv[])
 

Detailed Description

Analyse gravity wave data for tropical cyclones.

Definition in file hurricane.c.

Macro Definition Documentation

◆ NSTORM

#define NSTORM   9000

Definition at line 33 of file hurricane.c.

◆ NTIME

#define NTIME   140

Definition at line 36 of file hurricane.c.

Function Documentation

◆ get_storm_pos()

int get_storm_pos ( int  nobs,
double  time_wmo[NTIME],
double  lon_wmo[NTIME],
double  lat_wmo[NTIME],
double  wind_wmo[NTIME],
double  pres_wmo[NTIME],
double  t,
int  dt,
int  st,
double  x[3],
double *  wind,
double *  dwind,
double *  pres,
double *  dpres 
)

Definition at line 366 of file hurricane.c.

380 {
381
382 double x0[3], x1[3];
383
384 /* Check time range... */
385 if (t < time_wmo[0] || t > time_wmo[nobs - 1])
386 return 0;
387
388 /* Interpolate position... */
389 int i = locate_irr(time_wmo, nobs, t);
390 double w = (t - time_wmo[i]) / (time_wmo[i + 1] - time_wmo[i]);
391 geo2cart(0, lon_wmo[i], lat_wmo[i], x0);
392 geo2cart(0, lon_wmo[i + 1], lat_wmo[i + 1], x1);
393 x[0] = (1 - w) * x0[0] + w * x1[0];
394 x[1] = (1 - w) * x0[1] + w * x1[1];
395 x[2] = (1 - w) * x0[2] + w * x1[2];
396
397 /* Interpolate wind and pressure... */
398 *pres = (1 - w) * pres_wmo[i] + w * pres_wmo[i + 1];
399 *wind = (1 - w) * wind_wmo[i] + w * wind_wmo[i + 1];
400
401 /* Get pressure and wind change... */
402 *dpres = (pres_wmo[i + 1 + st] - pres_wmo[GSL_MAX(i - dt + st, 0)])
403 / (time_wmo[i + 1 + st] - time_wmo[GSL_MAX(i - dt + st, 0)]) * 3600.;
404 *dwind = (wind_wmo[i + 1 + st] - wind_wmo[GSL_MAX(i - dt + st, 0)])
405 / (time_wmo[i + 1 + st] - time_wmo[GSL_MAX(i - dt + st, 0)]) * 3600.;
406
407 return 1;
408}
locate_irr
int locate_irr(const double *xx, const int n, const double x)
Find array index for irregular grid.
Definition: jurassic.c:4103
geo2cart
void geo2cart(const double z, const double lon, const double lat, double *x)
Convert geolocation to Cartesian coordinates.
Definition: jurassic.c:3500
Here is the call graph for this function:

◆ read_var()

void read_var ( int  ncid,
const char  varname[],
size_t  nstorm,
int  nobs[NSTORM],
double  x[NSTORM][NTIME] 
)

Definition at line 412 of file hurricane.c.

417 {
418
419 int varid;
420
421 size_t count[2], start[2];
422
423 /* Read pressure... */
424 NC(nc_inq_varid(ncid, varname, &varid));
425 for (size_t istorm = 0; istorm < nstorm; istorm++) {
426 start[0] = istorm;
427 start[1] = 0;
428 count[0] = 1;
429 count[1] = (size_t) nobs[istorm];
430 NC(nc_get_vara_double(ncid, varid, start, count, x[istorm]));
431 }
432}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: diff_apr.c:35

◆ main()

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

Definition at line 71 of file hurricane.c.

73 {
74
75 static pert_t *pert;
76
77 static FILE *out;
78
79 static char filter[LEN], pertname[LEN], set[LEN];
80
81 static double bt4_mean, bt4_var, bt8_min, dpres, dpresbest, dwind,
82 dwindbest, lat_wmo[NSTORM][NTIME], latbest, lon_wmo[NSTORM][NTIME],
83 lonbest, lonsat, lonstorm, nedt, pres_wmo[NSTORM][NTIME],
84 pres, presbest, r2best = 1e100, wind_wmo[NSTORM][NTIME], wind,
85 windbest, time_max_pres[NSTORM], time_max_wind[NSTORM],
86 time_wmo[NSTORM][NTIME], timebest, xf[PERT_NTRACK][PERT_NXTRACK][3],
87 xs[3], z;
88
89 static int dimid, n, ncid, nobs[NSTORM], varid;
90
91 static size_t nstorm, ntime;
92
93 /* Check arguments... */
94 if (argc < 5)
95 ERRMSG("Give parameters: <ctl> <hurr.tab> <ibtracs.nc>"
96 " <pert1.nc> [<pert2.nc> ...]");
97
98 /* Get control parameters... */
99 scan_ctl(argc, argv, "SET", -1, "full", set);
100 scan_ctl(argc, argv, "PERTNAME", -1, "4mu", pertname);
101 scan_ctl(argc, argv, "FILTER", -1, "both", filter);
102 const double dt230 = scan_ctl(argc, argv, "DT230", -1, "0.16", NULL);
103 const double nu = scan_ctl(argc, argv, "NU", -1, "2345.0", NULL);
104 const double rmax = scan_ctl(argc, argv, "RMAX", -1, "500", NULL);
105 const int dt = (int) scan_ctl(argc, argv, "DT", -1, "0", NULL);
106 const int st = (int) scan_ctl(argc, argv, "ST", -1, "0", NULL);
107
108 /* Allocate... */
109 ALLOC(pert, pert_t, 1);
110
111 /* ------------------------------------------------------------
112 Read hurricane tracks...
113 ------------------------------------------------------------ */
114
115 /* Write info... */
116 printf("Read hurricane tracks: %s\n", argv[3]);
117
118 /* Open netCDF file... */
119 NC(nc_open(argv[3], NC_NOWRITE, &ncid));
120
121 /* Get dimensions... */
122 NC(nc_inq_dimid(ncid, "storm", &dimid));
123 NC(nc_inq_dimlen(ncid, dimid, &nstorm));
124 NC(nc_inq_dimid(ncid, "time", &dimid));
125 NC(nc_inq_dimlen(ncid, dimid, &ntime));
126 if (nstorm > NSTORM)
127 ERRMSG("Too many storms!");
128 if (ntime > NTIME)
129 ERRMSG("Too many time steps!");
130
131 /* Read number of observations per storm... */
132 NC(nc_inq_varid(ncid, "numObs", &varid));
133 NC(nc_get_var_int(ncid, varid, nobs));
134
135 /* Read data... */
136 read_var(ncid, "lat_wmo", nstorm, nobs, lat_wmo);
137 read_var(ncid, "lon_wmo", nstorm, nobs, lon_wmo);
138 read_var(ncid, "time_wmo", nstorm, nobs, time_wmo);
139 read_var(ncid, "wind_wmo", nstorm, nobs, wind_wmo);
140 read_var(ncid, "pres_wmo", nstorm, nobs, pres_wmo);
141
142 /* Convert units.. */
143 for (size_t istorm = 0; istorm < nstorm; istorm++)
144 for (int iobs = 0; iobs < nobs[istorm]; iobs++) {
145 time_wmo[istorm][iobs] *= 86400.;
146 time_wmo[istorm][iobs] -= 4453401600.00;
147 lon_wmo[istorm][iobs] *= 0.01;
148 lat_wmo[istorm][iobs] *= 0.01;
149 wind_wmo[istorm][iobs] *= 0.0514444;
150 pres_wmo[istorm][iobs] *= 0.1;
151 }
152
153 /* Check data... */
154 for (size_t istorm = 0; istorm < nstorm; istorm++)
155 for (int iobs = 0; iobs < nobs[istorm]; iobs++) {
156 if (pres_wmo[istorm][iobs] <= 800 || pres_wmo[istorm][iobs] >= 1200)
157 pres_wmo[istorm][iobs] = GSL_NAN;
158 if (wind_wmo[istorm][iobs] <= 0.1)
159 wind_wmo[istorm][iobs] = GSL_NAN;
160 }
161
162 /* Find time of maximum intensity (lowest pressure)... */
163 for (size_t istorm = 0; istorm < nstorm; istorm++) {
164 double pmin = 1e100;
165 time_max_pres[istorm] = GSL_NAN;
166 for (int iobs = 0; iobs < nobs[istorm]; iobs++)
167 if (gsl_finite(pres_wmo[istorm][iobs]) && pres_wmo[istorm][iobs] < pmin) {
168 pmin = pres_wmo[istorm][iobs];
169 time_max_pres[istorm] = time_wmo[istorm][iobs];
170 }
171 }
172
173 /* Find time of maximum intensity (maximum wind)... */
174 for (size_t istorm = 0; istorm < nstorm; istorm++) {
175 double wmax = -1e100;
176 time_max_wind[istorm] = GSL_NAN;
177 for (int iobs = 0; iobs < nobs[istorm]; iobs++)
178 if (gsl_finite(wind_wmo[istorm][iobs]) && wind_wmo[istorm][iobs] > wmax) {
179 wmax = wind_wmo[istorm][iobs];
180 time_max_wind[istorm] = time_wmo[istorm][iobs];
181 }
182 }
183
184 /* Close netCDF file... */
185 NC(nc_close(ncid));
186
187 /* ------------------------------------------------------------
188 Analyze AIRS data...
189 ------------------------------------------------------------ */
190
191 /* Create file... */
192 printf("Write hurricane data: %s\n", argv[2]);
193 if (!(out = fopen(argv[2], "w")))
194 ERRMSG("Cannot create file!");
195
196 /* Write header... */
197 fprintf(out,
198 "# $1 = storm number\n"
199 "# $2 = storm time since first report [hr]\n"
200 "# $3 = storm time since wind maximum [hr]\n"
201 "# $4 = storm time since pressure minimum [hr]\n"
202 "# $5 = match time [s]\n"
203 "# $6 = match longitude [deg]\n"
204 "# $7 = match latitude [deg]\n"
205 "# $8 = match distance [km]\n"
206 "# $9 = wind speed [m/s]\n"
207 "# $10 = wind speed change [m/s/hr]\n");
208 fprintf(out,
209 "# $11 = pressure [hPa]\n"
210 "# $12 = pressure change [hPa/hr]\n"
211 "# $13 = 8.1 micron BT minimum [K]\n"
212 "# $14 = 4.3 micron BT variance [K^2]\n"
213 "# $15 = 4.3 micron BT variance (noise-corrected) [K^2]\n"
214 "# $16 = number of footprints\n\n");
215
216 /* Loop over perturbation files... */
217 for (int iarg = 4; iarg < argc; iarg++) {
218
219 /* Read perturbation data... */
220 FILE *in;
221 if (!(in = fopen(argv[iarg], "r")))
222 continue;
223 else {
224 fclose(in);
225 read_pert(argv[iarg], pertname, pert);
226 }
227
228 /* Get Cartesian coordinates... */
229 for (int itrack2 = 0; itrack2 < pert->ntrack; itrack2++)
230 for (int ixtrack2 = 0; ixtrack2 < pert->nxtrack; ixtrack2++)
231 geo2cart(0, pert->lon[itrack2][ixtrack2],
232 pert->lat[itrack2][ixtrack2], xf[itrack2][ixtrack2]);
233
234 /* Loop over storms... */
235 for (size_t istorm = 0; istorm < nstorm; istorm++) {
236
237 /* Loop along AIRS center track... */
238 for (int itrack = 0; itrack < pert->ntrack; itrack++) {
239
240 /* Get storm position... */
241 if (get_storm_pos(nobs[istorm], time_wmo[istorm], lon_wmo[istorm],
242 lat_wmo[istorm], wind_wmo[istorm], pres_wmo[istorm],
243 pert->time[itrack][pert->nxtrack / 2], dt, st, xs,
244 &wind, &dwind, &pres, &dpres)) {
245
246 /* Get distance... */
247 double r2 = DIST2(xs, xf[itrack][pert->nxtrack / 2]);
248
249 /* Find best match... */
250 if (r2 < r2best) {
251
252 /* Save position... */
253 r2best = r2;
254 timebest = pert->time[itrack][pert->nxtrack / 2];
255 cart2geo(xs, &z, &lonbest, &latbest);
256
257 /* Save wind... */
258 windbest = wind;
259 dwindbest = dwind;
260 presbest = pres;
261 dpresbest = dpres;
262
263 /* Get BT data... */
264 n = 0;
265 bt8_min = 1e100;
266 bt4_mean = 0;
267 bt4_var = 0;
268 for (int itrack2 = GSL_MAX(itrack - ((int)(rmax / 17) + 1), 0);
269 itrack2 <= GSL_MIN(itrack + ((int) (rmax / 17) + 1),
270 pert->ntrack - 1); itrack2++)
271 for (int ixtrack2 = 0; ixtrack2 < pert->nxtrack; ixtrack2++) {
272
273 /* Check data... */
274 if (pert->time[itrack2][ixtrack2] < 0
275 || pert->lon[itrack2][ixtrack2] < -180
276 || pert->lon[itrack2][ixtrack2] > 180
277 || pert->lat[itrack2][ixtrack2] < -90
278 || pert->lat[itrack2][ixtrack2] > 90
279 || pert->pt[itrack2][ixtrack2] < -100
280 || pert->pt[itrack2][ixtrack2] > 100
281 || !gsl_finite(pert->bt[itrack2][ixtrack2])
282 || !gsl_finite(pert->pt[itrack2][ixtrack2])
283 || !gsl_finite(pert->var[itrack2][ixtrack2])
284 || !gsl_finite(pert->dc[itrack2][ixtrack2]))
285 continue;
286
287 /* Check east/west filter... */
288 lonsat = pert->lon[itrack2][ixtrack2];
289 while (lonsat < 20)
290 lonsat += 360;
291 lonstorm = lonbest;
292 while (lonstorm < 20)
293 lonstorm += 360;
294 if ((filter[0] == 'e' || filter[0] == 'E')
295 && lonsat < lonstorm)
296 continue;
297 if ((filter[0] == 'w' || filter[0] == 'W')
298 && lonsat > lonstorm)
299 continue;
300
301 /* Get distance... */
302 if (DIST2(xs, xf[itrack2][ixtrack2]) < rmax * rmax) {
303 bt8_min = GSL_MIN(bt8_min, pert->dc[itrack2][ixtrack2]);
304 bt4_mean += pert->bt[itrack2][ixtrack2];
305 bt4_var += gsl_pow_2(pert->pt[itrack2][ixtrack2]);
306 n++;
307 }
308 }
309 }
310 }
311
312 /* Output over poles... */
313 if (fabs(pert->lat[itrack][pert->nxtrack / 2]) > 80.) {
314
315 /* Get and check ascending/descending flag... */
316 const int asc =
317 (pert->lat[itrack > 0 ? itrack : itrack + 1][pert->nxtrack / 2]
318 > pert->lat[itrack >
319 0 ? itrack - 1 : itrack][pert->nxtrack / 2]);
320 if ((set[0] == 'f' || set[0] == 'F')
321 || ((set[0] == 'a' || set[0] == 'A') && asc)
322 || ((set[0] == 'd' || set[0] == 'D') && !asc)) {
323
324 /* Check for match... */
325 if (r2best < 890. * 890.) {
326
327 /* Estimate noise... */
328 if (dt230 > 0) {
329 const double nesr =
330 PLANCK(230.0 + dt230, nu) - PLANCK(230.0, nu);
331 nedt =
332 BRIGHT(PLANCK(bt4_mean / n, nu) + nesr, nu) - bt4_mean / n;
333 }
334
335 /* Write output... */
336 if (n > 0)
337 fprintf(out,
338 "%lu %g %g %g %.2f %g %g %g %g %g %g %g %g %g %g %d\n",
339 istorm, (timebest - time_wmo[istorm][0]) / 3600.,
340 (timebest - time_max_wind[istorm]) / 3600.,
341 (timebest - time_max_pres[istorm]) / 3600.,
342 timebest, lonbest, latbest, sqrt(r2best), windbest,
343 dwindbest, presbest, dpresbest, bt8_min, bt4_var / n,
344 bt4_var / n - gsl_pow_2(nedt), n);
345 }
346 }
347
348 /* Reset... */
349 r2best = 1e100;
350 }
351 }
352 }
353 }
354
355 /* Close file... */
356 fclose(out);
357
358 /* Free... */
359 free(pert);
360
361 return EXIT_SUCCESS;
362}
get_storm_pos
int get_storm_pos(int nobs, double time_wmo[NTIME], double lon_wmo[NTIME], double lat_wmo[NTIME], double wind_wmo[NTIME], double pres_wmo[NTIME], double t, int dt, int st, double x[3], double *wind, double *dwind, double *pres, double *dpres)
Definition: hurricane.c:366
NSTORM
#define NSTORM
Definition: hurricane.c:33
read_var
void read_var(int ncid, const char varname[], size_t nstorm, int nobs[NSTORM], double x[NSTORM][NTIME])
Definition: hurricane.c:412
NTIME
#define NTIME
Definition: hurricane.c:36
cart2geo
void cart2geo(const double *x, double *z, double *lon, double *lat)
Convert Cartesian coordinates to geolocation.
Definition: jurassic.c:108
scan_ctl
double scan_ctl(int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
Definition: jurassic.c:5114
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:383
BRIGHT
#define BRIGHT(rad, nu)
Compute brightness temperature.
Definition: jurassic.h:126
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory.
Definition: jurassic.h:121
PLANCK
#define PLANCK(T, nu)
Compute Planck function.
Definition: jurassic.h:183
DIST2
#define DIST2(a, b)
Compute squared distance between two vectors.
Definition: jurassic.h:137
read_pert
void read_pert(char *filename, char *pertname, pert_t *pert)
Read radiance perturbation data.
Definition: libairs.c:1137
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libairs.h:126
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libairs.h:123
pert_t
Perturbation data.
Definition: libairs.h:205
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libairs.h:232
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libairs.h:229
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libairs.h:214
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libairs.h:223
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libairs.h:208
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libairs.h:211
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK]
Longitude [deg].
Definition: libairs.h:217
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (4 or 15 micron) [K].
Definition: libairs.h:226
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK]
Latitude [deg].
Definition: libairs.h:220
Here is the call graph for this function:
Generated on Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4