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 42 of file invert.c.

44 {
45
46 static ctl_t ctl;
47
48 static atm_t atm, atm2;
49
50 static obs_t obs;
51
52 static FILE *in, *out;
53
54 static char line[LEN];
55
56 static double rtime[NLMAX], rz[NLMAX], rlon[NLMAX], rlat[NLMAX], obs_meas,
57 obs_sim, scl = 1.0, scl_err, c0, c1, cov00, cov01, cov11, sumsq,
58 x[NMAX], x2[NMAX], y[NMAX], y_err[NMAX], y2[NMAX], y2_err[NMAX],
59 y2_sim[NMAX], y2_sim_err[NMAX], w2[NMAX];
60
61 static float rp[NLMAX], rt[NLMAX], rso2[NLMAX], rh2o[NLMAX],
62 ro3[NLMAX], robs[NLMAX];
63
64 static int i, ig, n, nl, ndata[NMAX], nprof;
65
66 /* Check arguments... */
67 if (argc < 6)
68 ERRMSG("Give parameters: <ctl> <prof> <inv> <atm> <kernel>");
69
70 /* Read control parameters... */
71 read_ctl(argc, argv, &ctl);
72 const double dt = scan_ctl(argc, argv, "INVERT_DT", -1, "86400", NULL);
73 const double obs_err =
74 scan_ctl(argc, argv, "INVERT_OBS_ERR", -1, "1.0", NULL);
75 const int data = (int) scan_ctl(argc, argv, "INVERT_DATA", -1, "2", NULL);
76 const int fit = (int) scan_ctl(argc, argv, "INVERT_FIT", -1, "3", NULL);
77 const int itmax =
78 (int) scan_ctl(argc, argv, "INVERT_ITMAX", -1, "10", NULL);
79 const double tol = scan_ctl(argc, argv, "INVERT_TOL", -1, "1e-4", NULL);
80
81 /* Check control parameters... */
82 if (ctl.ng != 4)
83 ERRMSG("Set NG = 4!");
84 if (strcmp(ctl.emitter[0], "SO2") != 0)
85 ERRMSG("Set EMITTER[0] = SO2!");
86 if (strcmp(ctl.emitter[1], "H2O") != 0)
87 ERRMSG("Set EMITTER[1] = H2O!");
88 if (strcmp(ctl.emitter[2], "O3") != 0)
89 ERRMSG("Set EMITTER[2] = O3!");
90 if (strcmp(ctl.emitter[3], "CO2") != 0)
91 ERRMSG("Set EMITTER[3] = CO2!");
92 if (ctl.nd != 2)
93 ERRMSG("Set ND = 2!");
94
95 /* Set control parameters... */
96 ctl.write_bbt = 1;
97 ctl.write_matrix = 1;
98
99 /* Set observation data... */
100 obs.nr = 1;
101 obs.obsz[0] = 705;
102
103 /* Initialize look-up tables... */
104 tbl_t *tbl = read_tbl(&ctl);
105
106 /* ------------------------------------------------------------
107 Read profiles...
108 ------------------------------------------------------------ */
109
110 /* Read profile data... */
111 LOG(1, "Read profile data: %s", argv[2]);
112
113 /* Open file... */
114 if (!(in = fopen(argv[2], "r")))
115 ERRMSG("Cannot open file!");
116
117 /* Read file... */
118 while (fgets(line, LEN, in)) {
119
120 /* Read data... */
121 if (sscanf(line, "%lg %lg %lg %lg %g %g %g %g %g %g", &rtime[nl],
122 &rz[nl], &rlon[nl], &rlat[nl], &rp[nl], &rt[nl], &rso2[nl],
123 &rh2o[nl], &ro3[nl], &robs[nl]) != 10)
124 continue;
125 if ((++nl) > NLMAX)
126 ERRMSG("Too many profile data points!");
127 }
128
129 /* Close files... */
130 fclose(in);
131
132 /* ------------------------------------------------------------
133 Fit scaling factor for total mass...
134 ------------------------------------------------------------ */
135
136 /* Iterations... */
137 for (int it = 0; it < itmax; it++) {
138
139 /* Init... */
140 atm.np = n = 0;
141 for (i = 0; i < NMAX; i++) {
142 ndata[i] = 0;
143 x[i] = y[i] = NAN;
144 }
145
146 /* Loop over lines... */
147 for (int il = 0; il < nl; il++) {
148
149 /* Check for new profile... */
150 if ((rtime[il] != atm.time[0]
151 || rlon[il] != atm.lon[0]
152 || rlat[il] != atm.lat[0])
153 && atm.np > 0) {
154
155 /* Call forward model... */
156 formod(&ctl, tbl, &atm, &obs);
157 obs_sim = obs.rad[0][0] - obs.rad[1][0];
158
159 /* Get time index... */
160 i = (int) ((atm.time[0] - rtime[0]) / dt);
161 if (i < 0 && i >= NMAX)
162 ERRMSG("Time index out of range!");
163
164 /* Get maxima... */
165 if (data == 1) {
166 x[i] = (isfinite(x[i]) ? MAX(x[i], obs_sim) : obs_sim);
167 y[i] = (isfinite(y[i]) ? MAX(y[i], obs_meas) : obs_meas);
168 y_err[i] = obs_err;
169 if (isfinite(x[i]) && isfinite(y[i]))
170 ndata[i] = 1;
171 }
172
173 /* Get means... */
174 else if (data == 2) {
175 if (ndata[i] == 0) {
176 x[i] = obs_sim;
177 y[i] = obs_meas;
178 y_err[i] = POW2(obs_meas);
179 } else {
180 x[i] += obs_sim;
181 y[i] += obs_meas;
182 y_err[i] += POW2(obs_meas);
183 }
184 ndata[i]++;
185 }
186
187 /* Calculate mean atmospheric profile... */
188 nprof++;
189 atm2.np = atm.np;
190 for (int ip = 0; ip < atm.np; ip++) {
191 atm2.time[ip] += atm.time[ip];
192 atm2.z[ip] += atm.z[ip];
193 atm2.lon[ip] += atm.lon[ip];
194 atm2.lat[ip] += atm.lat[ip];
195 atm2.p[ip] += atm.p[ip];
196 atm2.t[ip] += atm.t[ip];
197 for (ig = 1; ig < ctl.ng; ig++)
198 atm2.q[ig][ip] += atm.q[ig][ip];
199 }
200
201 /* Reset counter... */
202 atm.np = 0;
203 }
204
205 /* Save data... */
206 obs_meas = robs[il];
207 atm.time[atm.np] = rtime[il];
208 atm.z[atm.np] = rz[il];
209 atm.lon[atm.np] = rlon[il];
210 atm.lat[atm.np] = rlat[il];
211 atm.p[atm.np] = rp[il];
212 atm.t[atm.np] = rt[il];
213 atm.q[0][atm.np] = rso2[il] * scl;
214 atm.q[1][atm.np] = rh2o[il];
215 atm.q[2][atm.np] = ro3[il];
216 atm.q[3][atm.np] = 371.789948e-6 + 2.026214e-6
217 * (atm.time[atm.np] - 63158400.) / 31557600.;
218 if ((++atm.np) > NP)
219 ERRMSG("Too many data points!");
220 }
221
222 /* Calculate means... */
223 if (data == 2)
224 for (i = 0; i < NMAX; i++)
225 if (ndata[i] > 0) {
226 x[i] /= ndata[i];
227 y[i] /= ndata[i];
228 y_err[i] = sqrt(MAX(y_err[i] / ndata[i] - POW2(y[i]), 0.0))
229 / sqrt(ndata[i]); /* standard error! */
230 }
231
232 /* Filter data... */
233 n = 0;
234 for (i = 0; i < NMAX; i++)
235 if (ndata[i] > 0 && isfinite(x[i]) && isfinite(y[i])
236 && isfinite(y_err[i])) {
237 x2[n] = x[i];
238 y2[n] = y[i];
239 y2_err[n] = y_err[i];
240 w2[n] = 1. / POW2(y_err[i]);
241 n++;
242 }
243
244 /* Fit radiance data... */
245 if (fit == 1)
246 gsl_fit_mul(x2, 1, y2, 1, (size_t) n, &c1, &cov11, &sumsq);
247 else if (fit == 2)
248 gsl_fit_wmul(x2, 1, w2, 1, y2, 1, (size_t) n, &c1, &cov11, &sumsq);
249 else if (fit == 3)
250 gsl_fit_linear(x2, 1, y2, 1, (size_t) n, &c0, &c1, &cov00, &cov01,
251 &cov11, &sumsq);
252 else if (fit == 4)
253 gsl_fit_wlinear(x2, 1, w2, 1, y2, 1, (size_t) n, &c0, &c1, &cov00,
254 &cov01, &cov11, &sumsq);
255 else
256 ERRMSG("Check INVERT_FIT!");
257
258 /* Get new scaling factor... */
259 double scl_old = scl;
260 scl_err = scl * sqrt(cov11);
261 scl *= c1;
262
263 /* Write info... */
264 LOG(1, " it= %d | scl= %g +/- %g | RMSE= %g",
265 it, scl, scl_err, sqrt(sumsq / n));
266
267 /* Convergence test... */
268 if (fabs(2.0 * (scl - scl_old) / (scl + scl_old)) < tol)
269 break;
270 }
271
272 /* ------------------------------------------------------------
273 Write inversion data...
274 ------------------------------------------------------------ */
275
276 /* Write info... */
277 LOG(1, "Write inversion data: %s", argv[3]);
278
279 /* Create file... */
280 if (!(out = fopen(argv[3], "w")))
281 ERRMSG("Cannot create file!");
282
283 /* Write header... */
284 fprintf(out,
285 "# $1 = time (seconds since 2000-01-01T00:00Z)\n"
286 "# $2 = simulated SO2 index [K]\n"
287 "# $3 = scaled simulated SO2 index [K]\n"
288 "# $4 = error of scaled simulated SO2 index [K]\n"
289 "# $5 = observed SO2 index [K]\n"
290 "# $6 = error of observed SO2 index [K]\n\n");
291
292 /* Write data... */
293 for (i = 0; i < n; i++) {
294
295 /* Calculate scaled SO2 index... */
296 if (fit == 1 || fit == 2)
297 gsl_fit_mul_est(x2[i], c1, cov11, &y2_sim[i], &y2_sim_err[i]);
298 else if (fit == 3 || fit == 4)
299 gsl_fit_linear_est(x2[i], c0, c1, cov00, cov01, cov11, &y2_sim[i],
300 &y2_sim_err[i]);
301
302 /* Write output... */
303 fprintf(out, "%.2f %g %g %g %g %g\n", rtime[0] + (i + 0.5) * dt,
304 x2[i], y2_sim[i], y2_sim_err[i], y2[i], y2_err[i]);
305 }
306
307 /* Report scaling factor for total mass... */
308 fprintf(out, "\n");
309 fprintf(out, "# scl= %g +/- %g\n", scl, scl_err);
310 fprintf(out, "# c1= %g +/- %g\n", c1, sqrt(cov11));
311 if (fit == 3 || fit == 4) {
312 fprintf(out, "# c0= %g +/- %g\n", c0, sqrt(cov00));
313 fprintf(out, "# corr= %g\n", cov01 / (sqrt(cov00) * sqrt(cov11)));
314 }
315 fprintf(out, "# RMSE= %g\n", sqrt(sumsq / n));
316 fprintf(out, "# n= %d\n", n);
317
318 /* Close files... */
319 fclose(out);
320
321 /* ------------------------------------------------------------
322 Calculate kernel...
323 ------------------------------------------------------------ */
324
325 /* Set atmospheric data... */
326 for (int ip = 0; ip < atm2.np; ip++) {
327 atm2.time[ip] /= nprof;
328 atm2.z[ip] /= nprof;
329 atm2.lon[ip] /= nprof;
330 atm2.lat[ip] /= nprof;
331 atm2.p[ip] /= nprof;
332 atm2.t[ip] /= nprof;
333 for (ig = 0; ig < ctl.ng; ig++)
334 atm2.q[ig][ip] /= nprof;
335 }
336
337 /* Get sizes... */
338 const size_t nk = atm2x(&ctl, &atm2, NULL, NULL, NULL);
339 const size_t mk = obs2y(&ctl, &obs, NULL, NULL, NULL);
340
341 /* Allocate... */
342 gsl_matrix *k = gsl_matrix_alloc(mk, nk);
343
344 /* Compute kernel matrix... */
345 kernel(&ctl, tbl, &atm2, &obs, k);
346
347 /* Write atmospheric data... */
348 write_atm(NULL, argv[4], &ctl, &atm);
349
350 /* Write matrix to file... */
351 write_matrix(NULL, argv[5], &ctl, k, &atm2, &obs, "y", "x", "r");
352
353 /* Free... */
354 gsl_matrix_free(k);
355 free(tbl);
356
357 return EXIT_SUCCESS;
358}
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)
Write atmospheric data.
Definition: jurassic.c:5470
read_ctl
void read_ctl(int argc, char *argv[], ctl_t *ctl)
Read forward model control parameters.
Definition: jurassic.c:4658
kernel
void kernel(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute Jacobians.
Definition: jurassic.c:4128
scan_ctl
double scan_ctl(int argc, char *argv[], const char *varname, const int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
Definition: jurassic.c:5248
formod
void formod(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs)
Determine ray paths and compute radiative transfer.
Definition: jurassic.c:3204
read_tbl
tbl_t * read_tbl(const ctl_t *ctl)
Read look-up table data.
Definition: jurassic.c:5074
obs2y
size_t obs2y(const ctl_t *ctl, const obs_t *obs, gsl_vector *y, int *ida, int *ira)
Compose measurement vector.
Definition: jurassic.c:4294
atm2x
size_t atm2x(const ctl_t *ctl, const atm_t *atm, gsl_vector *x, int *iqa, int *ipa)
Compose state vector or parameter vector.
Definition: jurassic.c:29
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)
Write matrix.
Definition: jurassic.c:5625
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:394
POW2
#define POW2(x)
Compute x^2.
Definition: jurassic.h:188
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:238
NP
#define NP
Maximum number of atmospheric data points.
Definition: jurassic.h:374
LOG
#define LOG(level,...)
Print log message.
Definition: jurassic.h:222
MAX
#define MAX(a, b)
Macro to determine the maximum of two values.
Definition: jurassic.h:157
atm_t
Atmospheric data.
Definition: jurassic.h:494
atm_t::time
double time[NP]
Time (seconds since 2000-01-01T00:00Z).
Definition: jurassic.h:500
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: jurassic.h:509
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: jurassic.h:506
atm_t::t
double t[NP]
Temperature [K].
Definition: jurassic.h:515
atm_t::np
int np
Number of data points.
Definition: jurassic.h:497
atm_t::z
double z[NP]
Altitude [km].
Definition: jurassic.h:503
atm_t::q
double q[NG][NP]
Volume mixing ratio [ppv].
Definition: jurassic.h:518
atm_t::p
double p[NP]
Pressure [hPa].
Definition: jurassic.h:512
ctl_t
Forward model control parameters.
Definition: jurassic.h:547
ctl_t::write_matrix
int write_matrix
Write matrix file (0=no, 1=yes).
Definition: jurassic.h:688
ctl_t::ng
int ng
Number of emitters.
Definition: jurassic.h:550
ctl_t::nd
int nd
Number of radiance channels.
Definition: jurassic.h:568
ctl_t::emitter
char emitter[NG][LEN]
Name of each emitter.
Definition: jurassic.h:553
ctl_t::write_bbt
int write_bbt
Use brightness temperature instead of radiance (0=no, 1=yes).
Definition: jurassic.h:685
obs_t
Observation geometry and radiance data.
Definition: jurassic.h:761
obs_t::rad
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
Definition: jurassic.h:800
obs_t::obsz
double obsz[NR]
Observer altitude [km].
Definition: jurassic.h:770
obs_t::nr
int nr
Number of ray paths.
Definition: jurassic.h:764
tbl_t
Emissivity look-up tables.
Definition: jurassic.h:805
Here is the call graph for this function:
Generated on Sun Nov 2 2025 18:29:33 for JURASSIC by doxygen 1.9.4