Perform forward model calculations in a single directory.
136 {
137
138 static atm_t atm, atm2;
139 static obs_t obs, obs2;
140
141
142 read_obs(wrkdir, obsfile, ctl, &obs);
143
144
145 read_atm(wrkdir, atmfile, ctl, &atm);
146
147
148 if (task[0] == 'p' || task[0] == 'P') {
149
150
151 for (
int ir = 0; ir < obs.
nr; ir++) {
152
153
155 for (
int ip = 0; ip < atm.
np; ip++)
158 atm2.
z[atm2.
np] = atm.
z[ip];
161 atm2.
p[atm2.
np] = atm.
p[ip];
162 atm2.
t[atm2.
np] = atm.
t[ip];
163 for (
int ig = 0; ig < ctl->
ng; ig++)
164 atm2.
q[ig][atm2.
np] = atm.
q[ig][ip];
165 for (
int iw = 0; iw < ctl->
nw; iw++)
166 atm2.
k[iw][atm2.
np] = atm.
k[iw][ip];
168 }
169
170
173 obs2.
vpz[0] = obs.
vpz[ir];
179
180
182
183
184 formod(ctl, tbl, &atm2, &obs2);
185
186
187 for (
int id = 0;
id < ctl->
nd;
id++) {
188 obs.
rad[id][ir] = obs2.
rad[id][0];
189 obs.
tau[id][ir] = obs2.
tau[id][0];
190 }
191 }
192 }
193
194
196 }
197
198
199 else {
200
201
202 formod(ctl, tbl, &atm, &obs);
203
204
206
207
208 if (task[0] == 'c' || task[0] == 'C') {
209
211
212
217
218
219 for (
int ig = 0; ig < ctl->
ng; ig++) {
220
221
223
224
225 for (
int iw = 0; iw < ctl->
nw; iw++)
226 for (
int ip = 0; ip < atm2.
np; ip++)
228
229
230 for (
int ig2 = 0; ig2 < ctl->
ng; ig2++)
231 if (ig2 != ig)
232 for (
int ip = 0; ip < atm2.
np; ip++)
234
235
236 formod(ctl, tbl, &atm2, &obs);
237
238
239 sprintf(filename,
"%s.%s", radfile, ctl->
emitter[ig]);
241 }
242
243
245
246
247 for (
int ig = 0; ig < ctl->
ng; ig++)
248 for (
int ip = 0; ip < atm2.
np; ip++)
250
251
252 formod(ctl, tbl, &atm2, &obs);
253
254
255 sprintf(filename, "%s.EXTINCT", radfile);
257 }
258
259
260 if (task[0] == 't' || task[0] == 'T') {
261
262
263 double t_min, t_max, t_mean = 0, t_sigma = 0;
264 int n = 0;
265
266
267 gsl_rng_env_setup();
268 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
269
270
271 do {
272
273
275 double dtemp = 40. * (gsl_rng_uniform(rng) - 0.5);
276 double dpress = 1. - 0.1 * gsl_rng_uniform(rng);
278 for (
int ig = 0; ig < ctl->
ng; ig++)
279 dq[ig] = 0.8 + 0.4 * gsl_rng_uniform(rng);
280 for (
int ip = 0; ip < atm2.
np; ip++) {
283 for (
int ig = 0; ig < ctl->
ng; ig++)
284 atm.
q[ig][ip] *= dq[ig];
285 }
286
287
288 double t0 = omp_get_wtime();
289 formod(ctl, tbl, &atm2, &obs);
290 double dt = omp_get_wtime() - t0;
291
292
293 t_mean += (dt);
295 if (n == 0 || dt < t_min)
296 t_min = dt;
297 if (n == 0 || dt > t_max)
298 t_max = dt;
299 n++;
300
301 } while (t_mean < 10.0);
302
303
304 t_mean /= (double) n;
305 t_sigma = sqrt(t_sigma / (
double) n -
POW2(t_mean));
306 printf("RUNTIME_MEAN = %g s\n", t_mean);
307 printf("RUNTIME_SIGMA = %g s\n", t_sigma);
308 printf("RUNTIME_MIN = %g s\n", t_min);
309 printf("RUNTIME_MAX = %g s\n", t_max);
310 printf(
"RAYS_PER_SECOND = %g", (
double) obs.
nr / t_mean);
311 }
312
313
314 if (task[0] == 's' || task[0] == 'S') {
315
316
319 formod(ctl, tbl, &atm, &obs);
321
322
323 for (double dz = 0.01; dz <= 2; dz *= 1.1) {
324 printf("STEPSIZE: \n");
325 for (double ds = 0.1; ds <= 50; ds *= 1.1) {
326
327
330
331
332 double t0 = omp_get_wtime();
333 formod(ctl, tbl, &atm, &obs);
334 double dt = omp_get_wtime() - t0;
335
336
337 double mean[
ND], sigma[
ND];
338 for (
int id = 0;
id < ctl->
nd;
id++) {
339 mean[id] = sigma[id] = 0;
340 int n = 0;
341 for (
int ir = 0; ir < obs.
nr; ir++) {
342 double err = 200. * (obs.
rad[id][ir] - obs2.
rad[id][ir])
343 / (obs.
rad[
id][ir] + obs2.
rad[
id][ir]);
344 mean[id] += err;
345 sigma[id] +=
POW2(err);
346 n++;
347 }
348 mean[id] /= n;
349 sigma[id] = sqrt(sigma[
id] / n -
POW2(mean[
id]));
350 }
351
352
353 printf("STEPSIZE: %g %g %g", ds, dz, dt);
354 for (
int id = 0;
id < ctl->
nd;
id++)
355 printf(" %g %g", mean[id], sigma[id]);
356 printf("\n");
357 }
358 }
359 }
360 }
361}
void read_atm(const char *dirname, const char *filename, const ctl_t *ctl, atm_t *atm)
Read atmospheric data.
void copy_obs(const ctl_t *ctl, obs_t *obs_dest, const obs_t *obs_src, const int init)
Copy and initialize observation data.
void read_obs(const char *dirname, const char *filename, const ctl_t *ctl, obs_t *obs)
Read observation data.
void write_obs(const char *dirname, const char *filename, const ctl_t *ctl, const obs_t *obs)
Write observation data.
void formod(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs)
Determine ray paths and compute radiative transfer.
void copy_atm(const ctl_t *ctl, atm_t *atm_dest, const atm_t *atm_src, const int init)
Copy and initialize atmospheric data.
#define LEN
Maximum length of ASCII data lines.
#define POW2(x)
Compute x^2.
#define ND
Maximum number of radiance channels.
#define NG
Maximum number of emitters.
double time[NP]
Time (seconds since 2000-01-01T00:00Z).
double k[NW][NP]
Extinction [km^-1].
double lat[NP]
Latitude [deg].
double lon[NP]
Longitude [deg].
double t[NP]
Temperature [K].
int np
Number of data points.
double z[NP]
Altitude [km].
double q[NG][NP]
Volume mixing ratio [ppv].
double p[NP]
Pressure [hPa].
int nw
Number of spectral windows.
int ctm_co2
Compute CO2 continuum (0=no, 1=yes).
int ctm_n2
Compute N2 continuum (0=no, 1=yes).
int ctm_h2o
Compute H2O continuum (0=no, 1=yes).
int ng
Number of emitters.
int nd
Number of radiance channels.
int ctm_o2
Compute O2 continuum (0=no, 1=yes).
char emitter[NG][LEN]
Name of each emitter.
double rayds
Maximum step length for raytracing [km].
double raydz
Vertical step length for raytracing [km].
Observation geometry and radiance data.
double tau[ND][NR]
Transmittance of ray path.
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
double vpz[NR]
View point altitude [km].
double vplat[NR]
View point latitude [deg].
double obslon[NR]
Observer longitude [deg].
double obslat[NR]
Observer latitude [deg].
double obsz[NR]
Observer altitude [km].
double vplon[NR]
View point longitude [deg].
double time[NR]
Time (seconds since 2000-01-01T00:00Z).
int nr
Number of ray paths.