JURASSIC
Functions
formod.c File Reference

JURASSIC forward model. More...

#include "jurassic.h"

Go to the source code of this file.

Functions

void call_formod (ctl_t *ctl, const tbl_t *tbl, const char *wrkdir, const char *obsfile, const char *atmfile, const char *radfile, const char *task)
 Perform forward model calculations in a single directory. More...
 
int main (int argc, char *argv[])
 

Detailed Description

JURASSIC forward model.

Definition in file formod.c.

Function Documentation

◆ call_formod()

void call_formod ( ctl_t ctl,
const tbl_t tbl,
const char *  wrkdir,
const char *  obsfile,
const char *  atmfile,
const char *  radfile,
const char *  task 
)

Perform forward model calculations in a single directory.

Definition at line 129 of file formod.c.

136 {
137
138 static atm_t atm, atm2;
139 static obs_t obs, obs2;
140
141 /* Read observation geometry... */
142 read_obs(wrkdir, obsfile, ctl, &obs);
143
144 /* Read atmospheric data... */
145 read_atm(wrkdir, atmfile, ctl, &atm);
146
147 /* Compute multiple profiles... */
148 if (task[0] == 'p' || task[0] == 'P') {
149
150 /* Loop over ray paths... */
151 for (int ir = 0; ir < obs.nr; ir++) {
152
153 /* Get atmospheric data... */
154 atm2.np = 0;
155 for (int ip = 0; ip < atm.np; ip++)
156 if (atm.time[ip] == obs.time[ir]) {
157 atm2.time[atm2.np] = atm.time[ip];
158 atm2.z[atm2.np] = atm.z[ip];
159 atm2.lon[atm2.np] = atm.lon[ip];
160 atm2.lat[atm2.np] = atm.lat[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];
167 atm2.np++;
168 }
169
170 /* Get observation data... */
171 obs2.nr = 1;
172 obs2.time[0] = obs.time[ir];
173 obs2.vpz[0] = obs.vpz[ir];
174 obs2.vplon[0] = obs.vplon[ir];
175 obs2.vplat[0] = obs.vplat[ir];
176 obs2.obsz[0] = obs.obsz[ir];
177 obs2.obslon[0] = obs.obslon[ir];
178 obs2.obslat[0] = obs.obslat[ir];
179
180 /* Check number of data points... */
181 if (atm2.np > 0) {
182
183 /* Call forward model... */
184 formod(ctl, tbl, &atm2, &obs2);
185
186 /* Save radiance data... */
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 /* Write radiance data... */
195 write_obs(wrkdir, radfile, ctl, &obs);
196 }
197
198 /* Compute single profile... */
199 else {
200
201 /* Call forward model... */
202 formod(ctl, tbl, &atm, &obs);
203
204 /* Save radiance data... */
205 write_obs(wrkdir, radfile, ctl, &obs);
206
207 /* Compute contributions... */
208 if (task[0] == 'c' || task[0] == 'C') {
209
210 char filename[LEN];
211
212 /* Switch off continua... */
213 ctl->ctm_co2 = 0;
214 ctl->ctm_h2o = 0;
215 ctl->ctm_n2 = 0;
216 ctl->ctm_o2 = 0;
217
218 /* Loop over emitters... */
219 for (int ig = 0; ig < ctl->ng; ig++) {
220
221 /* Copy atmospheric data... */
222 copy_atm(ctl, &atm2, &atm, 0);
223
224 /* Set extinction to zero... */
225 for (int iw = 0; iw < ctl->nw; iw++)
226 for (int ip = 0; ip < atm2.np; ip++)
227 atm2.k[iw][ip] = 0;
228
229 /* Set volume mixing ratios to zero... */
230 for (int ig2 = 0; ig2 < ctl->ng; ig2++)
231 if (ig2 != ig)
232 for (int ip = 0; ip < atm2.np; ip++)
233 atm2.q[ig2][ip] = 0;
234
235 /* Call forward model... */
236 formod(ctl, tbl, &atm2, &obs);
237
238 /* Save radiance data... */
239 sprintf(filename, "%s.%s", radfile, ctl->emitter[ig]);
240 write_obs(wrkdir, filename, ctl, &obs);
241 }
242
243 /* Copy atmospheric data... */
244 copy_atm(ctl, &atm2, &atm, 0);
245
246 /* Set volume mixing ratios to zero... */
247 for (int ig = 0; ig < ctl->ng; ig++)
248 for (int ip = 0; ip < atm2.np; ip++)
249 atm2.q[ig][ip] = 0;
250
251 /* Call forward model... */
252 formod(ctl, tbl, &atm2, &obs);
253
254 /* Save radiance data... */
255 sprintf(filename, "%s.EXTINCT", radfile);
256 write_obs(wrkdir, filename, ctl, &obs);
257 }
258
259 /* Measure CPU-time... */
260 if (task[0] == 't' || task[0] == 'T') {
261
262 /* Init... */
263 double t_min, t_max, t_mean = 0, t_sigma = 0;
264 int n = 0;
265
266 /* Initialize random number generator... */
267 gsl_rng_env_setup();
268 gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
269
270 /* Loop over profiles... */
271 do {
272
273 /* Create random atmosphere... */
274 copy_atm(ctl, &atm2, &atm, 0);
275 double dtemp = 40. * (gsl_rng_uniform(rng) - 0.5);
276 double dpress = 1. - 0.1 * gsl_rng_uniform(rng);
277 double dq[NG];
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++) {
281 atm.t[ip] += dtemp;
282 atm.p[ip] *= dpress;
283 for (int ig = 0; ig < ctl->ng; ig++)
284 atm.q[ig][ip] *= dq[ig];
285 }
286
287 /* Measure runtime... */
288 double t0 = omp_get_wtime();
289 formod(ctl, tbl, &atm2, &obs);
290 double dt = omp_get_wtime() - t0;
291
292 /* Get runtime statistics... */
293 t_mean += (dt);
294 t_sigma += POW2(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 /* Write results... */
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 /* Analyze effect of step size... */
314 if (task[0] == 's' || task[0] == 'S') {
315
316 /* Reference run... */
317 ctl->rayds = 0.1;
318 ctl->raydz = 0.01;
319 formod(ctl, tbl, &atm, &obs);
320 copy_obs(ctl, &obs2, &obs, 0);
321
322 /* Loop over step size... */
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 /* Set step size... */
328 ctl->rayds = ds;
329 ctl->raydz = dz;
330
331 /* Measure runtime... */
332 double t0 = omp_get_wtime();
333 formod(ctl, tbl, &atm, &obs);
334 double dt = omp_get_wtime() - t0;
335
336 /* Get differences... */
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 /* Write results... */
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.
Definition: jurassic.c:4581
void copy_obs(const ctl_t *ctl, obs_t *obs_dest, const obs_t *obs_src, const int init)
Copy and initialize observation data.
Definition: jurassic.c:3161
void read_obs(const char *dirname, const char *filename, const ctl_t *ctl, obs_t *obs)
Read observation data.
Definition: jurassic.c:4842
void write_obs(const char *dirname, const char *filename, const ctl_t *ctl, const obs_t *obs)
Write observation data.
Definition: jurassic.c:5834
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:3212
void copy_atm(const ctl_t *ctl, atm_t *atm_dest, const atm_t *atm_src, const int init)
Copy and initialize atmospheric data.
Definition: jurassic.c:3107
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:394
#define POW2(x)
Compute x^2.
Definition: jurassic.h:188
#define ND
Maximum number of radiance channels.
Definition: jurassic.h:364
#define NG
Maximum number of emitters.
Definition: jurassic.h:369
Atmospheric data.
Definition: jurassic.h:494
double time[NP]
Time (seconds since 2000-01-01T00:00Z).
Definition: jurassic.h:500
double k[NW][NP]
Extinction [km^-1].
Definition: jurassic.h:521
double lat[NP]
Latitude [deg].
Definition: jurassic.h:509
double lon[NP]
Longitude [deg].
Definition: jurassic.h:506
double t[NP]
Temperature [K].
Definition: jurassic.h:515
int np
Number of data points.
Definition: jurassic.h:497
double z[NP]
Altitude [km].
Definition: jurassic.h:503
double q[NG][NP]
Volume mixing ratio [ppv].
Definition: jurassic.h:518
double p[NP]
Pressure [hPa].
Definition: jurassic.h:512
int nw
Number of spectral windows.
Definition: jurassic.h:574
int ctm_co2
Compute CO2 continuum (0=no, 1=yes).
Definition: jurassic.h:607
int ctm_n2
Compute N2 continuum (0=no, 1=yes).
Definition: jurassic.h:613
int ctm_h2o
Compute H2O continuum (0=no, 1=yes).
Definition: jurassic.h:610
int ng
Number of emitters.
Definition: jurassic.h:550
int nd
Number of radiance channels.
Definition: jurassic.h:568
int ctm_o2
Compute O2 continuum (0=no, 1=yes).
Definition: jurassic.h:616
char emitter[NG][LEN]
Name of each emitter.
Definition: jurassic.h:553
double rayds
Maximum step length for raytracing [km].
Definition: jurassic.h:622
double raydz
Vertical step length for raytracing [km].
Definition: jurassic.h:625
Observation geometry and radiance data.
Definition: jurassic.h:761
double tau[ND][NR]
Transmittance of ray path.
Definition: jurassic.h:797
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
Definition: jurassic.h:800
double vpz[NR]
View point altitude [km].
Definition: jurassic.h:779
double vplat[NR]
View point latitude [deg].
Definition: jurassic.h:785
double obslon[NR]
Observer longitude [deg].
Definition: jurassic.h:773
double obslat[NR]
Observer latitude [deg].
Definition: jurassic.h:776
double obsz[NR]
Observer altitude [km].
Definition: jurassic.h:770
double vplon[NR]
View point longitude [deg].
Definition: jurassic.h:782
double time[NR]
Time (seconds since 2000-01-01T00:00Z).
Definition: jurassic.h:767
int nr
Number of ray paths.
Definition: jurassic.h:764
Here is the call graph for this function:

◆ main()

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

Definition at line 48 of file formod.c.

50 {
51
52 static ctl_t ctl;
53
54 /* Check arguments... */
55 if (argc < 5)
56 ERRMSG("Give parameters: <ctl> <obs> <atm> <rad>");
57
58 /* Read control parameters... */
59 read_ctl(argc, argv, &ctl);
60
61 /* Initialize look-up tables... */
62 tbl_t *tbl = read_tbl(&ctl);
63
64#ifdef UNIFIED
65
66 static atm_t atm;
67 static obs_t obs;
68
69 /* Read observation geometry... */
70 read_obs(NULL, argv[2], &ctl, &obs);
71
72 /* Read atmospheric data... */
73 read_atm(NULL, argv[3], &ctl, &atm);
74
75 /* Call forward model... */
76 jur_unified_init(argc, argv);
77 jur_unified_formod_multiple_packages(&atm, &obs, 1, NULL);
78
79 /* Save radiance data... */
80 write_obs(NULL, argv[4], &ctl, &obs);
81
82#else
83
84 char dirlist[LEN], task[LEN];
85
86 /* Get task... */
87 scan_ctl(argc, argv, "TASK", -1, "-", task);
88
89 /* Get dirlist... */
90 scan_ctl(argc, argv, "DIRLIST", -1, "-", dirlist);
91
92 /* Single forward calculation... */
93 if (dirlist[0] == '-')
94 call_formod(&ctl, tbl, NULL, argv[2], argv[3], argv[4], task);
95
96 /* Work on directory list... */
97 else {
98
99 /* Open directory list... */
100 FILE *in;
101 if (!(in = fopen(dirlist, "r")))
102 ERRMSG("Cannot open directory list!");
103
104 /* Loop over directories... */
105 char wrkdir[LEN];
106 while (fscanf(in, "%4999s", wrkdir) != EOF) {
107
108 /* Write info... */
109 LOG(1, "\nWorking directory: %s", wrkdir);
110
111 /* Call forward model... */
112 call_formod(&ctl, tbl, wrkdir, argv[2], argv[3], argv[4], task);
113 }
114
115 /* Close dirlist... */
116 fclose(in);
117 }
118
119#endif
120
121 /* Free... */
122 free(tbl);
123
124 return EXIT_SUCCESS;
125}
void call_formod(ctl_t *ctl, const tbl_t *tbl, const char *wrkdir, const char *obsfile, const char *atmfile, const char *radfile, const char *task)
Perform forward model calculations in a single directory.
Definition: formod.c:129
void read_ctl(int argc, char *argv[], ctl_t *ctl)
Read forward model control parameters.
Definition: jurassic.c:4686
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:5278
tbl_t * read_tbl(const ctl_t *ctl)
Read look-up table data.
Definition: jurassic.c:5104
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:238
#define LOG(level,...)
Print log message.
Definition: jurassic.h:222
Forward model control parameters.
Definition: jurassic.h:547
Emissivity look-up tables.
Definition: jurassic.h:805
Here is the call graph for this function: