JURASSIC
invert.c
Go to the documentation of this file.
1/*
2 This file is part of JURASSIC.
3
4 JURASSIC 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 JURASSIC 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 JURASSIC. If not, see <http://www.gnu.org/licenses/>.
16
17 Copyright (C) 2019-2026 Forschungszentrum Juelich GmbH
18*/
19
25#include "jurassic.h"
26#include <gsl/gsl_fit.h>
27
28/* ------------------------------------------------------------
29 Constants...
30 ------------------------------------------------------------ */
31
33#define NLMAX 30000000
34
36#define NMAX 1000
37
38/* ------------------------------------------------------------
39 Functions...
40 ------------------------------------------------------------ */
41
43static void usage(
44 void);
45
46/* ------------------------------------------------------------
47 Main...
48 ------------------------------------------------------------ */
49
50int main(
51 int argc,
52 char *argv[]) {
53
54 static ctl_t ctl;
55
56 static atm_t atm, atm2;
57
58 static obs_t obs;
59
60 static FILE *in, *out;
61
62 static char line[LEN];
63
64 static double rtime[NLMAX], rz[NLMAX], rlon[NLMAX], rlat[NLMAX], obs_meas,
65 obs_sim, scl = 1.0, scl_err, c0, c1, cov00, cov01, cov11, sumsq,
66 x[NMAX], x2[NMAX], y[NMAX], y_err[NMAX], y2[NMAX], y2_err[NMAX],
67 y2_sim[NMAX], y2_sim_err[NMAX], w2[NMAX];
68
69 static float rp[NLMAX], rt[NLMAX], rso2[NLMAX], rh2o[NLMAX],
70 ro3[NLMAX], robs[NLMAX];
71
72 static int i, ig, n, nl, ndata[NMAX], nprof;
73
74 /* Print usage information... */
75 USAGE;
76
77 /* Check arguments... */
78 if (argc < 6)
79 ERRMSG("Missing or invalid command-line arguments.\n\n"
80 "Usage: invert <ctl> <prof> <inv> <atm> <kernel> [KEY VALUE ...]\n\n"
81 "Use -h for full help.");
82
83 /* Read control parameters... */
84 read_ctl(argc, argv, &ctl);
85 const double dt = scan_ctl(argc, argv, "INVERT_DT", -1, "86400", NULL);
86 const double obs_err =
87 scan_ctl(argc, argv, "INVERT_OBS_ERR", -1, "1.0", NULL);
88 const int data = (int) scan_ctl(argc, argv, "INVERT_DATA", -1, "2", NULL);
89 const int fit = (int) scan_ctl(argc, argv, "INVERT_FIT", -1, "3", NULL);
90 const int itmax =
91 (int) scan_ctl(argc, argv, "INVERT_ITMAX", -1, "10", NULL);
92 const double tol = scan_ctl(argc, argv, "INVERT_TOL", -1, "1e-4", NULL);
93
94 /* Check control parameters... */
95 if (ctl.ng != 4)
96 ERRMSG("Set NG = 4!");
97 if (strcmp(ctl.emitter[0], "SO2") != 0)
98 ERRMSG("Set EMITTER[0] = SO2!");
99 if (strcmp(ctl.emitter[1], "H2O") != 0)
100 ERRMSG("Set EMITTER[1] = H2O!");
101 if (strcmp(ctl.emitter[2], "O3") != 0)
102 ERRMSG("Set EMITTER[2] = O3!");
103 if (strcmp(ctl.emitter[3], "CO2") != 0)
104 ERRMSG("Set EMITTER[3] = CO2!");
105 if (ctl.nd != 2)
106 ERRMSG("Set ND = 2!");
107
108 /* Set control parameters... */
109 ctl.write_bbt = 1;
110 ctl.write_matrix = 1;
111
112 /* Set observation data... */
113 obs.nr = 1;
114 obs.obsz[0] = 705;
115
116 /* Initialize look-up tables... */
117 tbl_t *tbl = read_tbl(&ctl);
118
119 /* ------------------------------------------------------------
120 Read profiles...
121 ------------------------------------------------------------ */
122
123 /* Read profile data... */
124 LOG(1, "Read profile data: %s", argv[2]);
125
126 /* Open file... */
127 if (!(in = fopen(argv[2], "r")))
128 ERRMSG("Cannot open file!");
129
130 /* Read file... */
131 while (fgets(line, LEN, in)) {
132
133 /* Read data... */
134 if (sscanf(line, "%lg %lg %lg %lg %g %g %g %g %g %g", &rtime[nl],
135 &rz[nl], &rlon[nl], &rlat[nl], &rp[nl], &rt[nl], &rso2[nl],
136 &rh2o[nl], &ro3[nl], &robs[nl]) != 10)
137 continue;
138 if ((++nl) > NLMAX)
139 ERRMSG("Too many profile data points!");
140 }
141
142 /* Close files... */
143 fclose(in);
144
145 /* ------------------------------------------------------------
146 Fit scaling factor for total mass...
147 ------------------------------------------------------------ */
148
149 /* Iterations... */
150 for (int it = 0; it < itmax; it++) {
151
152 /* Init... */
153 atm.np = n = 0;
154 for (i = 0; i < NMAX; i++) {
155 ndata[i] = 0;
156 x[i] = y[i] = NAN;
157 }
158
159 /* Loop over lines... */
160 for (int il = 0; il < nl; il++) {
161
162 /* Check for new profile... */
163 if ((rtime[il] != atm.time[0]
164 || rlon[il] != atm.lon[0]
165 || rlat[il] != atm.lat[0])
166 && atm.np > 0) {
167
168 /* Call forward model... */
169 formod(&ctl, tbl, &atm, &obs);
170 obs_sim = obs.rad[0][0] - obs.rad[1][0];
171
172 /* Get time index... */
173 i = (int) ((atm.time[0] - rtime[0]) / dt);
174 if (i < 0 && i >= NMAX)
175 ERRMSG("Time index out of range!");
176
177 /* Get maxima... */
178 if (data == 1) {
179 x[i] = (isfinite(x[i]) ? MAX(x[i], obs_sim) : obs_sim);
180 y[i] = (isfinite(y[i]) ? MAX(y[i], obs_meas) : obs_meas);
181 y_err[i] = obs_err;
182 if (isfinite(x[i]) && isfinite(y[i]))
183 ndata[i] = 1;
184 }
185
186 /* Get means... */
187 else if (data == 2) {
188 if (ndata[i] == 0) {
189 x[i] = obs_sim;
190 y[i] = obs_meas;
191 y_err[i] = POW2(obs_meas);
192 } else {
193 x[i] += obs_sim;
194 y[i] += obs_meas;
195 y_err[i] += POW2(obs_meas);
196 }
197 ndata[i]++;
198 }
199
200 /* Calculate mean atmospheric profile... */
201 nprof++;
202 atm2.np = atm.np;
203 for (int ip = 0; ip < atm.np; ip++) {
204 atm2.time[ip] += atm.time[ip];
205 atm2.z[ip] += atm.z[ip];
206 atm2.lon[ip] += atm.lon[ip];
207 atm2.lat[ip] += atm.lat[ip];
208 atm2.p[ip] += atm.p[ip];
209 atm2.t[ip] += atm.t[ip];
210 for (ig = 1; ig < ctl.ng; ig++)
211 atm2.q[ig][ip] += atm.q[ig][ip];
212 }
213
214 /* Reset counter... */
215 atm.np = 0;
216 }
217
218 /* Save data... */
219 obs_meas = robs[il];
220 atm.time[atm.np] = rtime[il];
221 atm.z[atm.np] = rz[il];
222 atm.lon[atm.np] = rlon[il];
223 atm.lat[atm.np] = rlat[il];
224 atm.p[atm.np] = rp[il];
225 atm.t[atm.np] = rt[il];
226 atm.q[0][atm.np] = rso2[il] * scl;
227 atm.q[1][atm.np] = rh2o[il];
228 atm.q[2][atm.np] = ro3[il];
229 atm.q[3][atm.np] = CO2_FIT(atm.time[atm.np]);
230 if ((++atm.np) > NP)
231 ERRMSG("Too many data points!");
232 }
233
234 /* Calculate means... */
235 if (data == 2)
236 for (i = 0; i < NMAX; i++)
237 if (ndata[i] > 0) {
238 x[i] /= ndata[i];
239 y[i] /= ndata[i];
240 y_err[i] = sqrt(MAX(y_err[i] / ndata[i] - POW2(y[i]), 0.0))
241 / sqrt(ndata[i]); /* standard error! */
242 }
243
244 /* Filter data... */
245 n = 0;
246 for (i = 0; i < NMAX; i++)
247 if (ndata[i] > 0 && isfinite(x[i]) && isfinite(y[i])
248 && isfinite(y_err[i])) {
249 x2[n] = x[i];
250 y2[n] = y[i];
251 y2_err[n] = y_err[i];
252 w2[n] = 1. / POW2(y_err[i]);
253 n++;
254 }
255
256 /* Fit radiance data... */
257 if (fit == 1)
258 gsl_fit_mul(x2, 1, y2, 1, (size_t) n, &c1, &cov11, &sumsq);
259 else if (fit == 2)
260 gsl_fit_wmul(x2, 1, w2, 1, y2, 1, (size_t) n, &c1, &cov11, &sumsq);
261 else if (fit == 3)
262 gsl_fit_linear(x2, 1, y2, 1, (size_t) n, &c0, &c1, &cov00, &cov01,
263 &cov11, &sumsq);
264 else if (fit == 4)
265 gsl_fit_wlinear(x2, 1, w2, 1, y2, 1, (size_t) n, &c0, &c1, &cov00,
266 &cov01, &cov11, &sumsq);
267 else
268 ERRMSG("Check INVERT_FIT!");
269
270 /* Get new scaling factor... */
271 double scl_old = scl;
272 scl_err = scl * sqrt(cov11);
273 scl *= c1;
274
275 /* Write info... */
276 LOG(1, " it= %d | scl= %g +/- %g | RMSE= %g",
277 it, scl, scl_err, sqrt(sumsq / n));
278
279 /* Convergence test... */
280 if (fabs(2.0 * (scl - scl_old) / (scl + scl_old)) < tol)
281 break;
282 }
283
284 /* ------------------------------------------------------------
285 Write inversion data...
286 ------------------------------------------------------------ */
287
288 /* Write info... */
289 LOG(1, "Write inversion data: %s", argv[3]);
290
291 /* Create file... */
292 if (!(out = fopen(argv[3], "w")))
293 ERRMSG("Cannot create file!");
294
295 /* Write header... */
296 fprintf(out,
297 "# $1 = time (seconds since 2000-01-01T00:00Z)\n"
298 "# $2 = simulated SO2 index [K]\n"
299 "# $3 = scaled simulated SO2 index [K]\n"
300 "# $4 = error of scaled simulated SO2 index [K]\n"
301 "# $5 = observed SO2 index [K]\n"
302 "# $6 = error of observed SO2 index [K]\n\n");
303
304 /* Write data... */
305 for (i = 0; i < n; i++) {
306
307 /* Calculate scaled SO2 index... */
308 if (fit == 1 || fit == 2)
309 gsl_fit_mul_est(x2[i], c1, cov11, &y2_sim[i], &y2_sim_err[i]);
310 else if (fit == 3 || fit == 4)
311 gsl_fit_linear_est(x2[i], c0, c1, cov00, cov01, cov11, &y2_sim[i],
312 &y2_sim_err[i]);
313
314 /* Write output... */
315 fprintf(out, "%.2f %g %g %g %g %g\n", rtime[0] + (i + 0.5) * dt,
316 x2[i], y2_sim[i], y2_sim_err[i], y2[i], y2_err[i]);
317 }
318
319 /* Report scaling factor for total mass... */
320 fprintf(out, "\n");
321 fprintf(out, "# scl= %g +/- %g\n", scl, scl_err);
322 fprintf(out, "# c1= %g +/- %g\n", c1, sqrt(cov11));
323 if (fit == 3 || fit == 4) {
324 fprintf(out, "# c0= %g +/- %g\n", c0, sqrt(cov00));
325 fprintf(out, "# corr= %g\n", cov01 / (sqrt(cov00) * sqrt(cov11)));
326 }
327 fprintf(out, "# RMSE= %g\n", sqrt(sumsq / n));
328 fprintf(out, "# n= %d\n", n);
329
330 /* Close files... */
331 fclose(out);
332
333 /* ------------------------------------------------------------
334 Calculate kernel...
335 ------------------------------------------------------------ */
336
337 /* Set atmospheric data... */
338 for (int ip = 0; ip < atm2.np; ip++) {
339 atm2.time[ip] /= nprof;
340 atm2.z[ip] /= nprof;
341 atm2.lon[ip] /= nprof;
342 atm2.lat[ip] /= nprof;
343 atm2.p[ip] /= nprof;
344 atm2.t[ip] /= nprof;
345 for (ig = 0; ig < ctl.ng; ig++)
346 atm2.q[ig][ip] /= nprof;
347 }
348
349 /* Get sizes... */
350 const size_t nk = atm2x(&ctl, &atm2, NULL, NULL, NULL);
351 const size_t mk = obs2y(&ctl, &obs, NULL, NULL, NULL);
352
353 /* Allocate... */
354 gsl_matrix *k = gsl_matrix_alloc(mk, nk);
355
356 /* Compute kernel matrix... */
357 kernel(&ctl, tbl, &atm2, &obs, k);
358
359 /* Write atmospheric data... */
360 write_atm(NULL, argv[4], &ctl, &atm, 0);
361
362 /* Write matrix to file... */
363 write_matrix(NULL, argv[5], &ctl, k, &atm2, &obs, "y", "x", "r", 0);
364
365 /* Free... */
366 gsl_matrix_free(k);
367 tbl_free(&ctl, tbl);
368
369 return EXIT_SUCCESS;
370}
371
372/*****************************************************************************/
373
374static void usage(
375 void) {
376 printf("\nJURASSIC inversion tool.\n\n");
377 printf("Run the MPTRAC inversion workflow using profile data, forward\n");
378 printf("model calculations, and kernel output.\n\n");
379 printf("Usage:\n");
380 printf(" invert <ctl> <prof> <inv> <atm> <kernel> [KEY VALUE ...]\n\n");
381 printf("Arguments:\n");
382 printf(" <ctl> Control file.\n");
383 printf(" <prof> Input profile data file.\n");
384 printf(" <inv> Output inversion summary file.\n");
385 printf(" <atm> Output atmospheric profile file.\n");
386 printf(" <kernel> Output kernel file.\n");
387 printf(" [KEY VALUE] Optional control parameters.\n\n");
388 printf("Tool-specific control parameters:\n");
389 printf(" INVERT_DT Time bin size.\n");
390 printf(" INVERT_OBS_ERR Observation error.\n");
391 printf(" INVERT_DATA Data reduction mode.\n");
392 printf(" INVERT_FIT Fit mode.\n");
393 printf(" INVERT_ITMAX Maximum number of iterations.\n");
394 printf(" INVERT_TOL Iteration tolerance.\n\n");
395 printf("Common control parameters:\n");
396 printf
397 (" TBLBASE, TBLFMT Lookup-table base name and format.\n");
398 printf(" ATMFMT, OBSFMT, MATRIXFMT Input/output file formats.\n");
399 printf(" NG, EMITTER[i] Active emitters.\n");
400 printf(" ND, NU[i], NW, WINDOW[i] Spectral channels and windows.\n");
401 printf
402 (" NCL, CLNU[i], NSF, SFNU[i] Cloud and surface spectral grids.\n");
403 printf(" RET*_ZMIN, RET*_ZMAX State-vector altitude limits.\n");
404 printf
405 (" WRITE_BBT, FORMOD Output units and forward-model selection.\n");
406 printf(" CTM_*, REFRAC Continua and refractivity.\n");
407 printf
408 (" RAYDS, RAYDZ, FOV Ray tracing and field of view.\n\n");
409 printf("Further information:\n");
410 printf(" Manual: https://slcs-jsc.github.io/jurassic/\n");
411}
usage
void usage(void)
Print command-line help.
Definition: formod.c:163
main
int main(int argc, char *argv[])
Definition: invert.c:50
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:7307
tbl_free
void tbl_free(const ctl_t *ctl, tbl_t *tbl)
Free lookup table and all internally allocated memory.
Definition: jurassic.c:7002
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:5516
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:7753
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:4387
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:6612
formod
void formod(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs)
Execute the selected forward model.
Definition: jurassic.c:3421
read_tbl
tbl_t * read_tbl(const ctl_t *ctl)
Read emissivity lookup tables from disk.
Definition: jurassic.c:6332
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:4627
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
jurassic.h
JURASSIC library declarations.
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:268
POW2
#define POW2(x)
Compute the square of a value.
Definition: jurassic.h:1018
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: jurassic.h:1325
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: jurassic.h:1206
CO2_FIT
#define CO2_FIT(TIME)
Climatological CO2 volume mixing ratio as a quadratic function of time.
Definition: jurassic.h:485
NP
#define NP
Maximum number of atmospheric data points.
Definition: jurassic.h:308
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: jurassic.h:1255
MAX
#define MAX(a, b)
Determine the maximum of two values.
Definition: jurassic.h:666
atm_t
Atmospheric profile data.
Definition: jurassic.h:1375
atm_t::time
double time[NP]
Time (seconds since 2000-01-01T00:00Z).
Definition: jurassic.h:1381
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: jurassic.h:1390
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: jurassic.h:1387
atm_t::t
double t[NP]
Temperature [K].
Definition: jurassic.h:1396
atm_t::np
int np
Number of data points.
Definition: jurassic.h:1378
atm_t::z
double z[NP]
Altitude [km].
Definition: jurassic.h:1384
atm_t::q
double q[NG][NP]
Volume mixing ratio [ppv].
Definition: jurassic.h:1399
atm_t::p
double p[NP]
Pressure [hPa].
Definition: jurassic.h:1393
ctl_t
Control parameters.
Definition: jurassic.h:1428
ctl_t::write_matrix
int write_matrix
Write matrix file (0=no, 1=yes).
Definition: jurassic.h:1572
ctl_t::ng
int ng
Number of emitters.
Definition: jurassic.h:1431
ctl_t::nd
int nd
Number of radiance channels.
Definition: jurassic.h:1449
ctl_t::emitter
char emitter[NG][LEN]
Name of each emitter.
Definition: jurassic.h:1434
ctl_t::write_bbt
int write_bbt
Use brightness temperature instead of radiance (0=no, 1=yes).
Definition: jurassic.h:1569
obs_t
Observation geometry and radiance data.
Definition: jurassic.h:1657
obs_t::rad
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
Definition: jurassic.h:1696
obs_t::obsz
double obsz[NR]
Observer altitude [km].
Definition: jurassic.h:1666
obs_t::nr
int nr
Number of ray paths.
Definition: jurassic.h:1660
tbl_t
Emissivity look-up tables.
Definition: jurassic.h:1842
Generated by doxygen 1.9.4