JURASSIC
Macros | Functions
invert.c File Reference

Inversion tool for MPTRAC. More...

#include "jurassic.h"
#include <gsl/gsl_fit.h>

Go to the source code of this file.

Macros

#define NLMAX   30000000
 Maximum number of data lines... More...
 
#define NMAX   1000
 Maximum number of data points... More...
 

Functions

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

Detailed Description

Inversion tool for MPTRAC.

Definition in file invert.c.

Macro Definition Documentation

◆ NLMAX

#define NLMAX   30000000

Maximum number of data lines...

Definition at line 33 of file invert.c.

◆ NMAX

#define NMAX   1000

Maximum number of data points...

Definition at line 36 of file invert.c.

Function Documentation

◆ main()

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

Definition at line 50 of file invert.c.

52 {
53
54 static ctl_t ctl;
55
56 static atm_t atm, atm2;
57
58 static obs_t obs, obs_scratch;
59 static los_t los_scratch;
60
61 static FILE *in, *out;
62
63 static char line[LEN];
64
65 static double rtime[NLMAX], rz[NLMAX], rlon[NLMAX], rlat[NLMAX], obs_meas,
66 scl = 1.0, scl_err, c0, c1, cov00, cov01, cov11, sumsq,
67 x[NMAX], x2[NMAX], y[NMAX], y_err[NMAX], y2[NMAX], y2_err[NMAX],
68 y2_sim[NMAX], y2_sim_err[NMAX], w2[NMAX];
69
70 static float rp[NLMAX], rt[NLMAX], rso2[NLMAX], rh2o[NLMAX],
71 ro3[NLMAX], robs[NLMAX];
72
73 static int i, ig, n, nl, ndata[NMAX], nprof;
74
75 /* Print usage information... */
76 USAGE;
77
78 /* Check arguments... */
79 if (argc < 6)
80 ERRMSG("Missing or invalid command-line arguments.\n\n"
81 "Usage: invert <ctl> <prof> <inv> <atm> <kernel> [KEY VALUE ...]\n\n"
82 "Use -h for full help.");
83
84 /* Read control parameters... */
85 read_ctl(argc, argv, &ctl);
86 const double dt = scan_ctl(argc, argv, "INVERT_DT", -1, "86400", NULL);
87 const double obs_err =
88 scan_ctl(argc, argv, "INVERT_OBS_ERR", -1, "1.0", NULL);
89 const int data = (int) scan_ctl(argc, argv, "INVERT_DATA", -1, "2", NULL);
90 const int fit = (int) scan_ctl(argc, argv, "INVERT_FIT", -1, "3", NULL);
91 const int itmax =
92 (int) scan_ctl(argc, argv, "INVERT_ITMAX", -1, "10", NULL);
93 const double tol = scan_ctl(argc, argv, "INVERT_TOL", -1, "1e-4", NULL);
94
95 /* Check control parameters... */
96 if (ctl.ng != 4)
97 ERRMSG("Set NG = 4!");
98 if (strcmp(ctl.emitter[0], "SO2") != 0)
99 ERRMSG("Set EMITTER[0] = SO2!");
100 if (strcmp(ctl.emitter[1], "H2O") != 0)
101 ERRMSG("Set EMITTER[1] = H2O!");
102 if (strcmp(ctl.emitter[2], "O3") != 0)
103 ERRMSG("Set EMITTER[2] = O3!");
104 if (strcmp(ctl.emitter[3], "CO2") != 0)
105 ERRMSG("Set EMITTER[3] = CO2!");
106 if (ctl.nd != 2)
107 ERRMSG("Set ND = 2!");
108
109 /* Set control parameters... */
110 ctl.write_bbt = 1;
111 ctl.write_matrix = 1;
112
113 /* Set observation data... */
114 obs.nr = 1;
115 obs.obsz[0] = 705;
116
117 /* Initialize look-up tables... */
118 tbl_t *tbl = read_tbl(&ctl);
119
120 /* ------------------------------------------------------------
121 Read profiles...
122 ------------------------------------------------------------ */
123
124 /* Read profile data... */
125 LOG(1, "Read profile data: %s", argv[2]);
126
127 /* Open file... */
128 if (!(in = fopen(argv[2], "r")))
129 ERRMSG("Cannot open file!");
130
131 /* Read file... */
132 while (fgets(line, LEN, in)) {
133
134 /* Read data... */
135 if (sscanf(line, "%lg %lg %lg %lg %g %g %g %g %g %g", &rtime[nl],
136 &rz[nl], &rlon[nl], &rlat[nl], &rp[nl], &rt[nl], &rso2[nl],
137 &rh2o[nl], &ro3[nl], &robs[nl]) != 10)
138 continue;
139 if ((++nl) > NLMAX)
140 ERRMSG("Too many profile data points!");
141 }
142
143 /* Close files... */
144 fclose(in);
145
146 /* ------------------------------------------------------------
147 Fit scaling factor for total mass...
148 ------------------------------------------------------------ */
149
150 /* Iterations... */
151 for (int it = 0; it < itmax; it++) {
152
153 /* Init... */
154 atm.np = n = 0;
155 for (i = 0; i < NMAX; i++) {
156 ndata[i] = 0;
157 x[i] = y[i] = NAN;
158 }
159
160 /* Loop over lines... */
161 for (int il = 0; il < nl; il++) {
162
163 /* Check for new profile... */
164 if ((rtime[il] != atm.time[0]
165 || rlon[il] != atm.lon[0]
166 || rlat[il] != atm.lat[0])
167 && atm.np > 0) {
168
169 /* Call forward model... */
170 if (formod(&ctl, tbl, &atm, &obs, &los_scratch, &obs_scratch)
171 != FORMOD_STATUS_OK)
172 ERRMSG("Forward model failed!");
173 const double obs_sim = obs.rad[0][0] - obs.rad[1][0];
174
175 /* Get time index... */
176 i = (int) ((atm.time[0] - rtime[0]) / dt);
177 if (i < 0 || i >= NMAX)
178 ERRMSG("Time index out of range!");
179
180 /* Get maxima... */
181 if (data == 1) {
182 x[i] = (isfinite(x[i]) ? MAX(x[i], obs_sim) : obs_sim);
183 y[i] = (isfinite(y[i]) ? MAX(y[i], obs_meas) : obs_meas);
184 y_err[i] = obs_err;
185 if (isfinite(x[i]) && isfinite(y[i]))
186 ndata[i] = 1;
187 }
188
189 /* Get means... */
190 else if (data == 2) {
191 if (ndata[i] == 0) {
192 x[i] = obs_sim;
193 y[i] = obs_meas;
194 y_err[i] = POW2(obs_meas);
195 } else {
196 x[i] += obs_sim;
197 y[i] += obs_meas;
198 y_err[i] += POW2(obs_meas);
199 }
200 ndata[i]++;
201 }
202
203 /* Calculate mean atmospheric profile... */
204 nprof++;
205 atm2.np = atm.np;
206 for (int ip = 0; ip < atm.np; ip++) {
207 atm2.time[ip] += atm.time[ip];
208 atm2.z[ip] += atm.z[ip];
209 atm2.lon[ip] += atm.lon[ip];
210 atm2.lat[ip] += atm.lat[ip];
211 atm2.p[ip] += atm.p[ip];
212 atm2.t[ip] += atm.t[ip];
213 for (ig = 1; ig < ctl.ng; ig++)
214 atm2.q[ig][ip] += atm.q[ig][ip];
215 }
216
217 /* Reset counter... */
218 atm.np = 0;
219 }
220
221 /* Save data... */
222 obs_meas = robs[il];
223 atm.time[atm.np] = rtime[il];
224 atm.z[atm.np] = rz[il];
225 atm.lon[atm.np] = rlon[il];
226 atm.lat[atm.np] = rlat[il];
227 atm.p[atm.np] = rp[il];
228 atm.t[atm.np] = rt[il];
229 atm.q[0][atm.np] = rso2[il] * scl;
230 atm.q[1][atm.np] = rh2o[il];
231 atm.q[2][atm.np] = ro3[il];
232 atm.q[3][atm.np] = CO2_FIT(atm.time[atm.np]);
233 if ((++atm.np) > NP)
234 ERRMSG("Too many data points!");
235 }
236
237 /* Calculate means... */
238 if (data == 2)
239 for (i = 0; i < NMAX; i++)
240 if (ndata[i] > 0) {
241 x[i] /= ndata[i];
242 y[i] /= ndata[i];
243 y_err[i] = sqrt(MAX(y_err[i] / ndata[i] - POW2(y[i]), 0.0))
244 / sqrt(ndata[i]); /* standard error! */
245 }
246
247 /* Filter data... */
248 n = 0;
249 for (i = 0; i < NMAX; i++)
250 if (ndata[i] > 0 && isfinite(x[i]) && isfinite(y[i])
251 && isfinite(y_err[i])) {
252 x2[n] = x[i];
253 y2[n] = y[i];
254 y2_err[n] = y_err[i];
255 w2[n] = 1. / POW2(y_err[i]);
256 n++;
257 }
258
259 /* Fit radiance data... */
260 if (fit == 1)
261 gsl_fit_mul(x2, 1, y2, 1, (size_t) n, &c1, &cov11, &sumsq);
262 else if (fit == 2)
263 gsl_fit_wmul(x2, 1, w2, 1, y2, 1, (size_t) n, &c1, &cov11, &sumsq);
264 else if (fit == 3)
265 gsl_fit_linear(x2, 1, y2, 1, (size_t) n, &c0, &c1, &cov00, &cov01,
266 &cov11, &sumsq);
267 else if (fit == 4)
268 gsl_fit_wlinear(x2, 1, w2, 1, y2, 1, (size_t) n, &c0, &c1, &cov00,
269 &cov01, &cov11, &sumsq);
270 else
271 ERRMSG("Check INVERT_FIT!");
272
273 /* Get new scaling factor... */
274 double scl_old = scl;
275 scl_err = scl * sqrt(cov11);
276 scl *= c1;
277
278 /* Write info... */
279 LOG(1, " it= %d | scl= %g +/- %g | RMSE= %g",
280 it, scl, scl_err, sqrt(sumsq / n));
281
282 /* Convergence test... */
283 if (fabs(2.0 * (scl - scl_old) / (scl + scl_old)) < tol)
284 break;
285 }
286
287 /* ------------------------------------------------------------
288 Write inversion data...
289 ------------------------------------------------------------ */
290
291 /* Write info... */
292 LOG(1, "Write inversion data: %s", argv[3]);
293
294 /* Create file... */
295 if (!(out = fopen(argv[3], "w")))
296 ERRMSG("Cannot create file!");
297
298 /* Write header... */
299 fprintf(out,
300 "# $1 = time (seconds since 2000-01-01T00:00Z)\n"
301 "# $2 = simulated SO2 index [K]\n"
302 "# $3 = scaled simulated SO2 index [K]\n"
303 "# $4 = error of scaled simulated SO2 index [K]\n"
304 "# $5 = observed SO2 index [K]\n"
305 "# $6 = error of observed SO2 index [K]\n\n");
306
307 /* Write data... */
308 for (i = 0; i < n; i++) {
309
310 /* Calculate scaled SO2 index... */
311 if (fit == 1 || fit == 2)
312 gsl_fit_mul_est(x2[i], c1, cov11, &y2_sim[i], &y2_sim_err[i]);
313 else if (fit == 3 || fit == 4)
314 gsl_fit_linear_est(x2[i], c0, c1, cov00, cov01, cov11, &y2_sim[i],
315 &y2_sim_err[i]);
316
317 /* Write output... */
318 fprintf(out, "%.2f %g %g %g %g %g\n", rtime[0] + (i + 0.5) * dt,
319 x2[i], y2_sim[i], y2_sim_err[i], y2[i], y2_err[i]);
320 }
321
322 /* Report scaling factor for total mass... */
323 fprintf(out, "\n");
324 fprintf(out, "# scl= %g +/- %g\n", scl, scl_err);
325 fprintf(out, "# c1= %g +/- %g\n", c1, sqrt(cov11));
326 if (fit == 3 || fit == 4) {
327 fprintf(out, "# c0= %g +/- %g\n", c0, sqrt(cov00));
328 fprintf(out, "# corr= %g\n", cov01 / (sqrt(cov00) * sqrt(cov11)));
329 }
330 fprintf(out, "# RMSE= %g\n", sqrt(sumsq / n));
331 fprintf(out, "# n= %d\n", n);
332
333 /* Close files... */
334 fclose(out);
335
336 /* ------------------------------------------------------------
337 Calculate kernel...
338 ------------------------------------------------------------ */
339
340 /* Set atmospheric data... */
341 for (int ip = 0; ip < atm2.np; ip++) {
342 atm2.time[ip] /= nprof;
343 atm2.z[ip] /= nprof;
344 atm2.lon[ip] /= nprof;
345 atm2.lat[ip] /= nprof;
346 atm2.p[ip] /= nprof;
347 atm2.t[ip] /= nprof;
348 for (ig = 0; ig < ctl.ng; ig++)
349 atm2.q[ig][ip] /= nprof;
350 }
351
352 /* Get sizes... */
353 const size_t nk = atm2x(&ctl, &atm2, NULL, NULL, NULL);
354 const size_t mk = obs2y(&ctl, &obs, NULL, NULL, NULL);
355
356 /* Allocate... */
357 gsl_matrix *k = gsl_matrix_alloc(mk, nk);
358
359 /* Compute kernel matrix... */
360 kernel(&ctl, tbl, &atm2, &obs, k);
361
362 /* Write atmospheric data... */
363 write_atm(NULL, argv[4], &ctl, &atm, 0);
364
365 /* Write matrix to file... */
366 write_matrix(NULL, argv[5], &ctl, k, &atm2, &obs, "y", "x", "r", 0);
367
368 /* Free... */
369 gsl_matrix_free(k);
370 tbl_free(&ctl, tbl);
371
372 return EXIT_SUCCESS;
373}
NMAX
#define NMAX
Maximum number of data points...
Definition: invert.c:36
NLMAX
#define NLMAX
Maximum number of data lines...
Definition: invert.c:33
write_atm
void write_atm(const char *dirname, const char *filename, const ctl_t *ctl, const atm_t *atm, int profile)
Write atmospheric data to a file.
Definition: jurassic.c:7308
tbl_free
void tbl_free(const ctl_t *ctl, tbl_t *tbl)
Free lookup table and all internally allocated memory.
Definition: jurassic.c:7003
read_ctl
void read_ctl(int argc, char *argv[], ctl_t *ctl)
Read model control parameters from command-line and configuration input.
Definition: jurassic.c:5542
write_matrix
void write_matrix(const char *dirname, const char *filename, const ctl_t *ctl, const gsl_matrix *matrix, const atm_t *atm, const obs_t *obs, const char *rowspace, const char *colspace, const char *sort, int dataset)
Write a fully annotated matrix (e.g., Jacobian or gain matrix) to file.
Definition: jurassic.c:7754
kernel
void kernel(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute the Jacobian (kernel) matrix by finite differences.
Definition: jurassic.c:4367
scan_ctl
double scan_ctl(int argc, char *argv[], const char *varname, const int arridx, const char *defvalue, char *value)
Scan control file or command-line arguments for a configuration variable.
Definition: jurassic.c:6613
read_tbl
tbl_t * read_tbl(const ctl_t *ctl)
Read emissivity lookup tables from disk.
Definition: jurassic.c:6333
obs2y
size_t obs2y(const ctl_t *ctl, const obs_t *obs, gsl_vector *y, int *ida, int *ira)
Convert observation radiances into a measurement vector.
Definition: jurassic.c:4635
atm2x
size_t atm2x(const ctl_t *ctl, const atm_t *atm, gsl_vector *x, int *iqa, int *ipa)
Convert atmospheric data to state vector elements.
Definition: jurassic.c:125
formod
int formod(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs, los_t *los_scratch, obs_t *obs_scratch)
Execute one complete forward-model calculation.
Definition: jurassic.c:3265
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:276
FORMOD_STATUS_OK
@ FORMOD_STATUS_OK
Definition: jurassic.h:159
POW2
#define POW2(x)
Compute the square of a value.
Definition: jurassic.h:1049
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: jurassic.h:1356
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: jurassic.h:1237
CO2_FIT
#define CO2_FIT(TIME)
Climatological CO2 volume mixing ratio as a quadratic function of time.
Definition: jurassic.h:493
NP
#define NP
Maximum number of atmospheric data points.
Definition: jurassic.h:316
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: jurassic.h:1286
MAX
#define MAX(a, b)
Determine the maximum of two values.
Definition: jurassic.h:674
atm_t
Atmospheric profile data.
Definition: jurassic.h:1406
atm_t::time
double time[NP]
Time (seconds since 2000-01-01T00:00Z).
Definition: jurassic.h:1412
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: jurassic.h:1421
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: jurassic.h:1418
atm_t::t
double t[NP]
Temperature [K].
Definition: jurassic.h:1427
atm_t::np
int np
Number of data points.
Definition: jurassic.h:1409
atm_t::z
double z[NP]
Altitude [km].
Definition: jurassic.h:1415
atm_t::q
double q[NG][NP]
Volume mixing ratio [ppv].
Definition: jurassic.h:1430
atm_t::p
double p[NP]
Pressure [hPa].
Definition: jurassic.h:1424
ctl_t
Control parameters.
Definition: jurassic.h:1459
ctl_t::write_matrix
int write_matrix
Write matrix file (0=no, 1=yes).
Definition: jurassic.h:1603
ctl_t::ng
int ng
Number of emitters.
Definition: jurassic.h:1462
ctl_t::nd
int nd
Number of radiance channels.
Definition: jurassic.h:1480
ctl_t::emitter
char emitter[NG][LEN]
Name of each emitter.
Definition: jurassic.h:1465
ctl_t::write_bbt
int write_bbt
Use brightness temperature instead of radiance (0=no, 1=yes).
Definition: jurassic.h:1600
los_t
Line-of-sight data.
Definition: jurassic.h:1626
obs_t
Observation geometry and radiance data.
Definition: jurassic.h:1688
obs_t::rad
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
Definition: jurassic.h:1727
obs_t::obsz
double obsz[NR]
Observer altitude [km].
Definition: jurassic.h:1697
obs_t::nr
int nr
Number of ray paths.
Definition: jurassic.h:1691
tbl_t
Emissivity look-up tables.
Definition: jurassic.h:1876
Here is the call graph for this function:
Generated by doxygen 1.9.4