NASA Logo
Ocean Color Science Software

ocssw V2022
gas_trans.c
Go to the documentation of this file.
1 /* =========================================================== */
2 /* Module gaseous_transmittance.c */
3 /* */
4 /* Computes sensor-specific transmittance through various */
5 /* atmospheric gases. */
6 /* */
7 /* B. Franz, NASA/OBPG, July 2006 */
8 /* */
9 /* adding the capability to read LUT with air mass factor */
10 /* as one additional dimension */
11 /* M. Zhang, SAIC/NASA, Nov. 2024 */
12 /* =========================================================== */
13 
14 #include "l12_proto.h"
15 #include <allocate3d.h>
16 #include "atrem_corl1.h"
17 
18 float get_wv_band_ratio(l1str *l1rec,int32_t ip,float window1, float absorp_band,float window2);
19 
20 int32_t model = 5;
21 static int amf;
23 static float *wvtbl = NULL;
24 static float *t_co2 = NULL;
25 static float *t_o2 = NULL;
26 static float *t_co = NULL;
27 static float *t_ch4 = NULL;
28 static float *t_n2o = NULL;
29 static float *cwv_all = NULL;
30 
31 static float *amf_mixed = NULL;
32 static float *amf_wv = NULL;
33 
34 static int index_amf_solz, index_amf_senz, index_amf_total;
35 static float ratio_solz, ratio_senz,ratio_total; //interpolation ratio for amf at direction of solz,senz and two-way
36 static float amf_solz,amf_senz,amf_total;
37 
38 void load_gas_tables(l1str *l1rec) {
39  /*
40  netcdf file that has water vapor transmittance as a function of wavelength and
41  cwv (6 profiles x number of bands x 220 water vapor value)
42  filename = /OCDATAROOT/sensor[/subsensor]/<sensorName>_gas_trans.nc
43  */
44 
45  /* This will be the netCDF ID for the file and data variable. */
46  int32_t ncid, varid;
47  int32_t num_water_vapors_id, num_models_id, num_wavelengths_id, num_airmass_id;
48  amf = 0;
49  num_airmass=1;
50 
51  /* Open the file */
52  if ((nc_open(input->gas_transmittance_file, NC_NOWRITE, &ncid)) != NC_NOERR) {
53  printf("-E- %s: Failed to open %s\n", __FILE__, input->gas_transmittance_file);
54  exit(EXIT_FAILURE);
55  }
56  if ((nc_inq_dimid(ncid, "n_air_mass_factor", &num_airmass_id)) == NC_NOERR) {
57  amf = 1;
58  if ((nc_inq_dimlen(ncid, num_airmass_id, &num_airmass)) != NC_NOERR) {
59  printf("-E- %s: Failed to read dimension n_water_vapor\n", __FILE__);
60  exit(EXIT_FAILURE);
61  }
62  }
63  if ((nc_inq_dimid(ncid, "n_water_vapor", &num_water_vapors_id)) == NC_NOERR) {
64  if ((nc_inq_dimlen(ncid, num_water_vapors_id, &num_water_vapors)) != NC_NOERR) {
65  printf("-E- %s: Failed to read dimension n_water_vapor\n", __FILE__);
66  exit(EXIT_FAILURE);
67  }
68  } else{
69  printf("-E- %s: Failed to find dimension n_water_vapor\n", __FILE__);
70  exit(EXIT_FAILURE);
71  }
72  if((nc_inq_dimid(ncid, "nmodels", &num_models_id)) == NC_NOERR){
73  if((nc_inq_dimlen(ncid, num_models_id, &num_models)) != NC_NOERR){
74  printf("-E- %s: Failed to read dimension nmodels\n", __FILE__);
75  exit(EXIT_FAILURE);
76  }
77  } else{
78  printf("-E- %s: Failed to find dimension nmodels\n", __FILE__);
79  exit(EXIT_FAILURE);
80  }
81 
82  if((nc_inq_dimid(ncid, "nwavelengths", &num_wavelengths_id)) == NC_NOERR){
83  if((nc_inq_dimlen(ncid, num_wavelengths_id, &num_wavelengths)) != NC_NOERR){
84  printf("-E- %s: Failed to read dimension nwavelengths\n", __FILE__);
85  exit(EXIT_FAILURE);
86  }
87  } else {
88  printf("-E- %s: Failed to find dimension nwavelengths\n", __FILE__);
89  exit(EXIT_FAILURE);
90  }
91 
92  /* Read the water vapor transmittance */
93  if ((wvtbl = (float *)malloc(num_models*num_wavelengths*num_airmass*num_water_vapors*sizeof(float))) == NULL) {
94  printf("-E- %s line %d : error allocating memory for water vapor transmittance table.\n",
95  __FILE__, __LINE__);
96  exit(1);
97  }
98  if ((nc_inq_varid(ncid, "water_vapor_transmittance", &varid)) == NC_NOERR) {
99  if ((nc_get_var_float(ncid, varid,wvtbl)) != NC_NOERR){
100  printf("-E- %s: failed to read water_vapor_transmittance from %s\n", __FILE__,
101  input->gas_transmittance_file);
102  exit(EXIT_FAILURE);
103  }
104  } else {
105  printf("-E- %s: '%s' does not have water_vapor_transmittance.\n",
106  __FILE__, input->gas_transmittance_file);
107  exit(EXIT_FAILURE);
108  }
109 
110  if (amf) {
111  if ((amf_mixed = (float *)malloc(num_airmass * sizeof(float))) == NULL) {
112  printf("Error: allocating memory for air mass factor mixed\n");
113  exit(EXIT_FAILURE);
114  }
115  if ((nc_inq_varid(ncid, "air_mass_factor_mixed", &varid)) == NC_NOERR) {
116  if ((nc_get_var_float(ncid, varid, amf_mixed)) != NC_NOERR) {
117  printf("-E- %s: failed to read air mass factor mixed from %s\n", __FILE__,
118  input->gas_transmittance_file);
119  exit(EXIT_FAILURE);
120  }
121  } else {
122  printf("-E- %s: '%s' does not have air mass factor mixed\n", __FILE__,
123  input->gas_transmittance_file);
124  exit(EXIT_FAILURE);
125  }
126 
127  if ((amf_wv = (float *)malloc(num_airmass * sizeof(float))) == NULL) {
128  printf("Error: allocating memory for air mass factor for water vapor\n");
129  exit(EXIT_FAILURE);
130  }
131  if ((nc_inq_varid(ncid, "air_mass_factor_wv", &varid)) == NC_NOERR) {
132  if ((nc_get_var_float(ncid, varid, amf_wv)) != NC_NOERR) {
133  printf("-E- %s: failed to read air mass factor for water vapor from %s\n", __FILE__,
134  input->gas_transmittance_file);
135  exit(EXIT_FAILURE);
136  }
137  } else {
138  printf("-E- %s: '%s' does not have air mass factor for water vapor\n", __FILE__,
139  input->gas_transmittance_file);
140  exit(EXIT_FAILURE);
141  }
142  }
143 
144  /* Read the water vapor table */
145  if ((cwv_all = (float *)malloc(num_water_vapors*sizeof(float))) == NULL) {
146  printf("Error: allocating memory for water vapor table\n");
147  exit(EXIT_FAILURE);
148  }
149  if ((nc_inq_varid(ncid, "water_vapor", &varid)) == NC_NOERR) {
150  if ((nc_get_var_float(ncid, varid, cwv_all)) != NC_NOERR) {
151  printf("-E- %s: failed to read water_vapor from %s\n", __FILE__, input->gas_transmittance_file);
152  exit(EXIT_FAILURE);
153  }
154  } else {
155  printf("-E- %s: '%s' does not have water_vapor\n",
156  __FILE__, input->gas_transmittance_file);
157  exit(EXIT_FAILURE);
158  }
159 
160  /* Read the carbon monoxide transmittance */
161  if ((t_co = (float *)malloc(num_wavelengths*num_airmass*sizeof(float))) == NULL) {
162  printf("-E- %s line %d : error allocating memory for carbon monoxide transmittance table.\n",
163  __FILE__, __LINE__);
164  exit(1);
165  }
166  if ((nc_inq_varid(ncid, "carbon_monoxide_transmittance", &varid)) == NC_NOERR) {
167  if ((nc_get_var_float(ncid, varid, t_co)) != NC_NOERR) {
168  printf("-E- %s: failed to read carbon_monoxide_transmittance from %s\n", __FILE__,
169  input->gas_transmittance_file);
170  exit(EXIT_FAILURE);
171  }
172  } else {
173  printf("-E- %s: '%s' does not have carbon_monoxide_transmittance\n",
174  __FILE__, input->gas_transmittance_file);
175  exit(EXIT_FAILURE);
176  }
177 
178  /* Read the carbon dioxide transmittance */
179  if ((t_co2 = (float *)malloc(num_wavelengths*num_airmass*sizeof(float))) == NULL) {
180  printf("-E- %s line %d : error allocating memory for carbon dioxide transmittance table.\n", __FILE__,
181  __LINE__);
182  exit(1);
183  }
184  if ((nc_inq_varid(ncid, "carbon_dioxide_transmittance", &varid)) == NC_NOERR) {
185  if ((nc_get_var_float(ncid, varid, t_co2)) != NC_NOERR) {
186  printf("-E- %s: failed to read carbon_dioxide_transmittance from %s\n", __FILE__,
187  input->gas_transmittance_file);
188  exit(EXIT_FAILURE);
189  }
190  } else {
191  printf("-E- %s: '%s' does not have carbon_dioxide_transmittance\n",
192  __FILE__, input->gas_transmittance_file);
193  exit(EXIT_FAILURE);
194  }
195 
196  /* Read the oxygen transmittance */
197  if ((t_o2 = (float *)malloc(num_wavelengths*num_airmass*sizeof(float))) == NULL) {
198  printf("-E- %s line %d : error allocating memory for oxygen transmittance table.\n", __FILE__,
199  __LINE__);
200  exit(1);
201  }
202  if ((nc_inq_varid(ncid, "oxygen_transmittance", &varid)) == NC_NOERR) {
203  if ((nc_get_var_float(ncid, varid, t_o2)) != NC_NOERR) {
204  printf("-E- %s: failed to read oxygen_transmittance from %s\n", __FILE__,
205  input->gas_transmittance_file);
206  exit(EXIT_FAILURE);
207  }
208  } else {
209  printf("-E- %s: '%s' does not have oxygen_transmittance\n",
210  __FILE__, input->gas_transmittance_file);
211  exit(EXIT_FAILURE);
212  }
213 
214  /* Read the nitrous oxide transmittance */
215  if ((t_n2o = (float *)malloc(num_wavelengths*num_airmass*sizeof(float))) == NULL) {
216  printf("-E- %s line %d : error allocating memory for nitrous oxide transmittance table.\n", __FILE__,
217  __LINE__);
218  exit(1);
219  }
220  if ((nc_inq_varid(ncid, "nitrous_oxide_transmittance", &varid)) == NC_NOERR) {
221  if ((nc_get_var_float(ncid, varid, t_n2o)) != NC_NOERR) {
222  printf("-E- %s: failed to read nitrous_oxide_transmittance from %s\n", __FILE__,
223  input->gas_transmittance_file);
224  exit(EXIT_FAILURE);
225  }
226  } else {
227  printf("-E- %s: '%s' does not have nitrous_oxide_transmittance\n",
228  __FILE__, input->gas_transmittance_file);
229  exit(EXIT_FAILURE);
230  }
231 
232  /* Read the methane transmittance */
233  if ((t_ch4 = (float *)malloc(num_wavelengths*num_airmass*sizeof(float))) == NULL) {
234  printf("-E- %s line %d : error allocating memory for methane transmittance table.\n", __FILE__,
235  __LINE__);
236  exit(1);
237  }
238  if ((nc_inq_varid(ncid, "methane_transmittance", &varid)) == NC_NOERR) {
239  if ((nc_get_var_float(ncid, varid, t_ch4)) != NC_NOERR) {
240  printf("-E- %s: failed to read methane_transmittance from %s\n", __FILE__,
241  input->gas_transmittance_file);
242  exit(EXIT_FAILURE);
243  }
244  } else {
245  printf("-E- %s: '%s' does not have methane_transmittance\n",
246  __FILE__, input->gas_transmittance_file);
247  exit(EXIT_FAILURE);
248  }
249 
250  /* Read the Nitrogen dioxide */
251  if ((nc_inq_varid(ncid, "k_no2", &varid)) == NC_NOERR) {
252  if ((nc_get_var_float(ncid, varid, l1rec->l1file->k_no2)) != NC_NOERR) {
253  printf("-E- %s: failed to read k_no2 from %s\n", __FILE__,
254  input->gas_transmittance_file);
255  exit(EXIT_FAILURE);
256  }
257  }
258  /* Read the ozone cross-section */
259  if ((nc_inq_varid(ncid, "k_oz", &varid)) == NC_NOERR) {
260  if ((nc_get_var_float(ncid, varid, l1rec->l1file->k_oz)) != NC_NOERR) {
261  printf("-E- %s: failed to read k_oz from %s\n", __FILE__,
262  input->gas_transmittance_file);
263  exit(EXIT_FAILURE);
264  }
265  }
266 
267  /* Close the file */
268  if ((nc_close(ncid)) != NC_NOERR){
269  printf("-E- %s: failed to close %s\n",
270  __FILE__, input->gas_transmittance_file);
271  exit(EXIT_FAILURE);
272  }
273 }
274 
275 int32_t get_index_lowerbound(float *table_val, int32_t num_val, float val) {
276  int32_t index;
277  int32_t i;
278 
279  for (i = 0; i < num_val; i++)
280  if (val < table_val[i])
281  break;
282  index=MAX(i-1,0);
283  index=MIN(i-1,num_val-2);
284  return index;
285 }
286 
287 
288 void ozone_transmittance(l1str *l1rec, int32_t ip) {
289  float tau_oz;
290  int32_t iw;
291 
292  filehandle *l1file = l1rec->l1file;
293  int32_t nwave = l1file->nbands;
294  int32_t ipb = ip*nwave;
295 
296  for (iw = 0; iw < nwave; iw++) {
297  tau_oz = l1rec->oz[ip] * l1file->k_oz[iw];
298  l1rec->tg_sol[ipb + iw] *= exp(-(tau_oz / l1rec->csolz[ip]));
299  l1rec->tg_sen[ipb + iw] *= exp(-(tau_oz / l1rec->csenz[ip]));
300  l1rec->tg[ipb + iw]*=exp(-tau_oz* (1.0/ l1rec->csolz[ip]+1.0/ l1rec->csenz[ip]));
301  }
302 }
303 
304 void co2_transmittance(l1str *l1rec, int32_t ip) {
305  int32_t iw;
306 
307  filehandle *l1file = l1rec->l1file;
308  int32_t nwave = l1file->nbands;
309  int32_t ipb = ip*nwave;
310 
311  if ((input->gas_opt & GAS_TRANS_TBL_BIT) == 0) {
312  if (t_co2 == NULL) {
313  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "t_co2", (void **) &t_co2);
314  }
315  }
316 
317  for (iw = 0; iw < nwave; iw++) {
318 
319  if(amf){
320  int32_t index=iw*num_airmass;
321  float t_co2_interp;
322 
323  t_co2_interp=t_co2[index+index_amf_solz]*(1-ratio_solz)+t_co2[index+index_amf_solz+1]*ratio_solz;
324  l1rec->tg_sol[ipb + iw] *= t_co2_interp;
325 
326  t_co2_interp=t_co2[index+index_amf_senz]*(1-ratio_senz)+t_co2[index+index_amf_senz+1]*ratio_senz;
327  l1rec->tg_sen[ipb + iw] *= t_co2_interp;
328 
329  t_co2_interp=t_co2[index+index_amf_total]*(1-ratio_total)+t_co2[index+index_amf_total+1]*ratio_total;
330  l1rec->tg[ipb + iw] *= t_co2_interp;
331  } else {
332  l1rec->tg_sol[ipb + iw] *= pow(t_co2[iw], amf_solz);
333  l1rec->tg_sen[ipb + iw] *= pow(t_co2[iw], amf_senz);
334  l1rec->tg [ipb + iw] *= pow(t_co2[iw], amf_total);
335  }
336  }
337 
338 }
339 
340 void co_transmittance(l1str *l1rec, int32_t ip) {
341  int32_t iw;
342 
343  filehandle *l1file = l1rec->l1file;
344  int32_t nwave = l1file->nbands;
345  int32_t ipb = ip*nwave;
346 
347  /* Only compute if gas transmittance table was requested */
348  if ((input->gas_opt & GAS_TRANS_TBL_BIT) != 0) {
349  for (iw = 0; iw < nwave; iw++) {
350  if (amf) {
351 
352  int32_t index=iw*num_airmass;
353 
354  float t_co_interp=t_co[index+index_amf_solz]*(1-ratio_solz)+t_co[index+index_amf_solz+1]*ratio_solz;
355  l1rec->tg_sol[ipb + iw] *= t_co_interp;
356 
357  t_co_interp=t_co[index+index_amf_senz]*(1-ratio_senz)+t_co[index+index_amf_senz+1]*ratio_senz;
358  l1rec->tg_sen[ipb + iw] *= t_co_interp;
359 
360  t_co_interp=t_co[index+index_amf_total]*(1-ratio_total)+t_co[index+index_amf_total+1]*ratio_total;
361  l1rec->tg[ipb + iw] *= t_co_interp;
362  } else {
363  l1rec->tg_sol[ipb + iw] *= pow(t_co[iw], amf_solz);
364  l1rec->tg_sen[ipb + iw] *= pow(t_co[iw], amf_senz);
365  l1rec->tg [ipb + iw] *= pow(t_co[iw], amf_total);
366  }
367  }
368  } else {
369  printf("-E- carbon monoxide transmittance can only be computed if gas_opt includes bit %d\n", GAS_TRANS_TBL_BIT);
370  exit(EXIT_FAILURE);
371 
372  }
373 }
374 
375 void ch4_transmittance(l1str *l1rec, int32_t ip) {
376  int32_t iw;
377 
378  filehandle *l1file = l1rec->l1file;
379  int32_t nwave = l1file->nbands;
380  int32_t ipb = ip*nwave;
381 
382  /* Only compute if gas transmittance table was requested */
383  if ((input->gas_opt & GAS_TRANS_TBL_BIT) != 0) {
384  for (iw = 0; iw < nwave; iw++) {
385  if (amf) {
386  int32_t index=iw*num_airmass;
387 
388  float t_ch4_interp=t_ch4[index+index_amf_solz]*(1-ratio_solz)+t_ch4[index+index_amf_solz+1]*ratio_solz;
389  l1rec->tg_sol[ipb + iw] *= t_ch4_interp;
390 
391  t_ch4_interp=t_ch4[index+index_amf_senz]*(1-ratio_senz)+t_ch4[index+index_amf_senz+1]*ratio_senz;
392  l1rec->tg_sen[ipb + iw] *= t_ch4_interp;
393 
394  t_ch4_interp=t_ch4[index+index_amf_total]*(1-ratio_total)+t_ch4[index+index_amf_total+1]*ratio_total;
395  l1rec->tg[ipb + iw] *= t_ch4_interp;
396  } else {
397  l1rec->tg_sol[ipb + iw] *= pow(t_ch4[iw], amf_solz);
398  l1rec->tg_sen[ipb + iw] *= pow(t_ch4[iw], amf_senz);
399  l1rec->tg [ipb + iw] *= pow(t_ch4[iw], amf_total);
400  }
401  }
402  } else {
403  printf("-E- methane transmittance can only be computed if gas_opt includes bit %d\n", GAS_TRANS_TBL_BIT);
404  exit(EXIT_FAILURE);
405  }
406 }
407 
408 void o2_transmittance(l1str *l1rec, int32_t ip) {
409  int32_t iw;
410 
411  filehandle *l1file = l1rec->l1file;
412  int32_t nwave = l1file->nbands;
413  int32_t ipb = ip*nwave;
414 
415  /* Only compute if gas transmittance table was requested */
416  if ((input->gas_opt & GAS_TRANS_TBL_BIT) != 0) {
417  for (iw = 0; iw < nwave; iw++) {
418  if (amf) {
419  int32_t index=iw*num_airmass;
420 
421  float t_o2_interp=t_o2[index+index_amf_solz]*(1-ratio_solz)+t_o2[index+index_amf_solz+1]*ratio_solz;
422  l1rec->tg_sol[ipb + iw] *= t_o2_interp;
423 
424  t_o2_interp=t_o2[index+index_amf_senz]*(1-ratio_senz)+t_o2[index+index_amf_senz+1]*ratio_senz;
425  l1rec->tg_sen[ipb + iw] *= t_o2_interp;
426 
427  t_o2_interp=t_o2[index+index_amf_total]*(1-ratio_total)+t_o2[index+index_amf_total+1]*ratio_total;
428  l1rec->tg[ipb + iw] *=t_o2_interp;
429  } else {
430  l1rec->tg_sol[ipb + iw] *= pow(t_o2[iw], amf_solz);
431  l1rec->tg_sen[ipb + iw] *= pow(t_o2[iw], amf_senz);
432  l1rec->tg [ipb + iw] *= pow(t_o2[iw], amf_total);
433  }
434  }
435  }
436 }
437 
438 void n2o_transmittance(l1str *l1rec, int32_t ip) {
439  int32_t iw;
440 
441  filehandle *l1file = l1rec->l1file;
442  int32_t nwave = l1file->nbands;
443  int32_t ipb = ip*nwave;
444 
445  /* Only compute if gas transmittance table was requested */
446  if ((input->gas_opt & GAS_TRANS_TBL_BIT) != 0) {
447  for (iw = 0; iw < nwave; iw++) {
448  if (amf) {
449  int32_t index=iw*num_airmass;
450 
451  float t_n2o_interp=t_n2o[index+index_amf_solz]*(1-ratio_solz)+t_n2o[index+index_amf_solz+1]*ratio_solz;
452  l1rec->tg_sol[ipb + iw] *= t_n2o_interp;
453 
454  t_n2o_interp=t_n2o[index+index_amf_senz]*(1-ratio_senz)+t_n2o[index+index_amf_senz+1]*ratio_senz;
455  l1rec->tg_sen[ipb + iw] *= t_n2o_interp;
456 
457  t_n2o_interp=t_n2o[index+index_amf_total]*(1-ratio_total)+t_n2o[index+index_amf_total+1]*ratio_total;
458  l1rec->tg[ipb + iw] *= t_n2o_interp;
459  } else {
460  l1rec->tg_sol[ipb + iw] *= pow(t_n2o[iw], amf_solz);
461  l1rec->tg_sen[ipb + iw] *= pow(t_n2o[iw], amf_senz);
462  l1rec->tg [ipb + iw] *= pow(t_n2o[iw], amf_total);
463  }
464  }
465  } else {
466  printf("-E- nitrous oxide transmittance can only be computed if gas_opt includes bit %d\n", GAS_TRANS_TBL_BIT);
467  exit(EXIT_FAILURE);
468  }
469 }
470 
471 void no2_transmittance(l1str *l1rec, int32_t ip) {
472 
473  float a_285, a_225;
474  float tau_to200;
475  float no2_tr200;
476  int32_t iw;
477 
478  filehandle *l1file = l1rec->l1file;
479  int32_t nwave = l1file->nbands;
480  int32_t ipb = ip*nwave;
481 
482  float sec0 = 1.0 / l1rec->csolz[ip];
483  float sec = 1.0 / l1rec->csenz[ip];
484 
485  if (l1rec->no2_tropo[ip] > 0.0)
486  /* compute tropo no2 above 200m (Z.Ahmad)
487  no2_tr200 = exp(12.6615 + 0.61676*log(no2_tropo));
488  new, location-dependent method */
489  no2_tr200 = l1rec->no2_frac[ip] * l1rec->no2_tropo[ip];
490  else
491  no2_tr200 = 0.0;
492 
493 
494  for (iw = 0; iw < nwave; iw++) {
495 
496  if (l1file->k_no2[iw] > 0.0) {
497 
498  a_285 = l1file->k_no2[iw] * (1.0 - 0.003 * (285.0 - 294.0));
499  a_225 = l1file->k_no2[iw] * (1.0 - 0.003 * (225.0 - 294.0));
500 
501  tau_to200 = a_285 * no2_tr200 + a_225 * l1rec->no2_strat[ip];
502 
503  l1rec->tg_sol[ipb + iw] *= exp(-(tau_to200 * sec0));
504  l1rec->tg_sen[ipb + iw] *= exp(-(tau_to200 * sec));
505  l1rec->tg [ipb + iw] *= exp(-(tau_to200 * (sec+sec0)));
506  }
507  }
508 }
509 
510 void h2o_transmittance(l1str *l1rec, int32_t ip) {
511  static float *a_h2o = NULL;
512  static float *b_h2o = NULL;
513  static float *c_h2o = NULL;
514  static float *d_h2o = NULL;
515  static float *e_h2o = NULL;
516  static float *f_h2o = NULL;
517  static float *g_h2o = NULL;
518  int32_t index;
519 
520  float t_h2o;
521  int32_t iw;
522 
523  filehandle *l1file = l1rec->l1file;
524  int32_t nwave = l1file->nbands;
525  int32_t ipb = ip*nwave;
526  float wv = l1rec->wv[ip];
527 
528  if (amf && input->watervapor_bands) {
529  wv = 0;
530  for (iw = 0; iw < input->nbands_watervapor;) {
531  wv += get_wv_band_ratio(l1rec, ip, input->watervapor_bands[iw], input->watervapor_bands[iw + 1],
532  input->watervapor_bands[iw + 2]);
533  iw += 3;
534  }
535  wv /= (input->nbands_watervapor / 3);
536  }
537 
538  // Apply water vapor transmittance only for the MSE and multi-band AC from the netcdf
539  // if (input->aer_opt == AERRHMSEPS || input->aer_opt == AERRHMSEPS_lin || input->aer_opt == AERRHSM) {
540  if ((input->gas_opt & GAS_TRANS_TBL_BIT) != 0) {
541  int32_t ja = 0,ja_sen=0,ja_sol=0,ja_total=0;
542  float ratio_wv;
543  float f00,f11,f01,f10;
544  float ratio_amf_senz,ratio_amf_solz,ratio_amf_total,tempratio;
545  int index_amf_wv_senz,index_amf_wv_solz,index_amf_wv_total;
546 
547  if (amf) {
548  index_amf_wv_senz = get_index_lowerbound(amf_wv, num_airmass, amf_senz);
549  index_amf_wv_solz = get_index_lowerbound(amf_wv, num_airmass, amf_solz);
550  index_amf_wv_total = get_index_lowerbound(amf_wv, num_airmass, amf_total);
551 
552  ratio_amf_senz = (amf_senz - amf_wv[index_amf_wv_senz]) /
553  (amf_wv[index_amf_wv_senz + 1] - amf_wv[index_amf_wv_senz]);
554  ratio_amf_solz = (amf_solz - amf_wv[index_amf_wv_solz]) /
555  (amf_wv[index_amf_wv_solz + 1] - amf_wv[index_amf_wv_solz]);
556  ratio_amf_total = (amf_total - amf_wv[index_amf_wv_total]) /
557  (amf_wv[index_amf_wv_total + 1] - amf_wv[index_amf_wv_total]);
558  }
559  ja = get_index_lowerbound(cwv_all, num_water_vapors, wv );
560  ja_sen = get_index_lowerbound(cwv_all, num_water_vapors, wv*amf_senz );
561  ja_sol = get_index_lowerbound(cwv_all, num_water_vapors, wv*amf_solz );
562  ja_total = get_index_lowerbound(cwv_all, num_water_vapors, wv*amf_total );
563 
564  ratio_wv=(wv -cwv_all[ja])/(cwv_all[ja+1]-cwv_all[ja]);
565 
566 
567  for (iw = 0; iw < nwave; iw++) {
568  if (amf) {
570 
571  f00 = wvtbl[index+index_amf_wv_solz*num_water_vapors+ja];
572  f10 = wvtbl[index+(index_amf_wv_solz+1)*num_water_vapors+ja];
573  f01 = wvtbl[index+index_amf_wv_solz*num_water_vapors+ja+1];
574  f11 = wvtbl[index+(index_amf_wv_solz+1)*num_water_vapors+ja+1];
575 
576  t_h2o = (1. - ratio_amf_solz)*(1. - ratio_wv) * f00 + ratio_amf_solz * ratio_wv * f11 + ratio_amf_solz * (1. - ratio_wv) * f10 + ratio_wv * (1. - ratio_amf_solz) * f01;
577  l1rec->tg_sol[ipb + iw] *= t_h2o;
578 
579  f00 = wvtbl[index+index_amf_wv_senz*num_water_vapors+ja];
580  f10 = wvtbl[index+(index_amf_wv_senz+1)*num_water_vapors+ja];
581  f01 = wvtbl[index+index_amf_wv_senz*num_water_vapors+ja+1];
582  f11 = wvtbl[index+(index_amf_wv_senz+1)*num_water_vapors+ja+1];
583 
584  t_h2o = (1. - ratio_amf_senz)*(1. - ratio_wv) * f00 + ratio_amf_senz * ratio_wv * f11 + ratio_amf_senz * (1. - ratio_wv) * f10 + ratio_wv * (1. - ratio_amf_senz) * f01;
585  l1rec->tg_sen[ipb + iw] *= t_h2o;
586 
587  f00 = wvtbl[index+index_amf_wv_total*num_water_vapors+ja];
588  f10 = wvtbl[index+(index_amf_wv_total+1)*num_water_vapors+ja];
589  f01 = wvtbl[index+index_amf_wv_total*num_water_vapors+ja+1];
590  f11 = wvtbl[index+(index_amf_wv_total+1)*num_water_vapors+ja+1];
591 
592  t_h2o = (1. - ratio_amf_total)*(1. - ratio_wv) * f00 + ratio_amf_total * ratio_wv * f11 + ratio_amf_total * (1. - ratio_wv) * f10 + ratio_wv * (1. - ratio_amf_total) * f01;
593  l1rec->tg[ipb + iw] *= t_h2o;
594  } else {
596 
597  tempratio=(wv*amf_solz -cwv_all[ja_sol])/(cwv_all[ja_sol+1]-cwv_all[ja_sol]);
598  t_h2o=wvtbl[index+ja_sol]*(1-tempratio)+wvtbl[index+ja_sol+1]*tempratio;
599  l1rec->tg_sol[ipb + iw] *= t_h2o;
600 
601  tempratio=(wv*amf_senz -cwv_all[ja_sen])/(cwv_all[ja_sen+1]-cwv_all[ja_sen]);
602  t_h2o=wvtbl[index+ja_sen]*(1-tempratio)+wvtbl[index+ja_sen+1]*tempratio;
603  l1rec->tg_sen[ipb + iw] *= t_h2o;
604 
605  tempratio=(wv*amf_total -cwv_all[ja_total])/(cwv_all[ja_total+1]-cwv_all[ja_total]);
606  t_h2o=wvtbl[index+ja_total]*(1-tempratio)+wvtbl[index+ja_total+1]*tempratio;
607  l1rec->tg[ipb + iw] *= t_h2o;
608  }
609  }
610  } // otherwise apply Zia's tabel from Bo-cai
611  else {
612  if (a_h2o == NULL) {
613  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "a_h2o", (void **) &a_h2o);
614  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "b_h2o", (void **) &b_h2o);
615  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "c_h2o", (void **) &c_h2o);
616  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "d_h2o", (void **) &d_h2o);
617  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "e_h2o", (void **) &e_h2o);
618  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "f_h2o", (void **) &f_h2o);
619  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "g_h2o", (void **) &g_h2o);
620  }
621 
622  for (iw = 0; iw < nwave; iw++) {
623  t_h2o = a_h2o[iw] + wv * (b_h2o[iw] + wv * (c_h2o[iw] + wv * (d_h2o[iw]
624  + wv * (e_h2o[iw] + wv * (f_h2o[iw] + wv * g_h2o[iw])))));
625  l1rec->tg_sol[ipb + iw] *= pow(t_h2o, 1.0 / l1rec->csolz[ip]);
626  l1rec->tg_sen[ipb + iw] *= pow(t_h2o, 1.0 / l1rec->csenz[ip]);
627  l1rec->tg [ipb + iw] *= pow(t_h2o, 1.0 / l1rec->csenz[ip]+1.0 / l1rec->csolz[ip]);
628  }
629  }
630  return;
631 }
632 
633 void gaseous_transmittance(l1str *l1rec, int32_t ip) {
634  static int32_t firstRun = TRUE;
635  int ib, ipb;
636  int nwave = l1rec->l1file->nbands;
637 
638  if ((input->gas_opt & ATREM_BIT) != 0) {
639  if (input->oxaband_opt == 1 && ((input->atrem_opt & ATREM_O2) != 0)){
640  printf("ATREM O2 correction is incompatible with the Ding and Gordon approach\n");
641  printf("Either unset the atrem_opt bit %d or set oxaband_opt=0\n",ATREM_O2);
642  exit(EXIT_FAILURE);
643  }
644  static float *rhot, *tg_tot;
645  float airmass, A;
646 
647  if (firstRun) {
648  if ((rhot = (float *) calloc(nwave, sizeof (float))) == NULL) {
649  printf("-E- : Error allocating memory to rhot\n");
650  exit(EXIT_FAILURE);
651  }
652  if ((tg_tot = (float *) calloc(nwave, sizeof (float))) == NULL) {
653  printf("-E- : Error allocating memory to tg_tot\n");
654  exit(EXIT_FAILURE);
655  }
656  firstRun = FALSE;
657  }
658  /* Copy surface reflectance to temp. var for atrem calculation */
659  for (ib = 0; ib < nwave; ib++) {
660  ipb = ip * nwave + ib;
661  rhot [ib] = M_PI * l1rec->Lt[ipb] / l1rec->Fo[ib] / l1rec->csolz[ip];
662  }
663 
664  get_atrem_cor(l1rec, ip, rhot, tg_tot, &l1rec->tg_sol[ip * nwave], &l1rec->tg_sen[ip * nwave]);
665  if (input->atrem_splitpaths == 0) {
666  airmass = 1.0 / l1rec->csolz[ip] + 1.0 / l1rec->csenz[ip];
667 
668  for (ib = 0; ib < nwave; ib++) {
669  ipb = ip * nwave + ib;
670  // ATREM didn't do the path splitting, do it here
671  A = -log(tg_tot[ib]) / airmass; //effective optical depth
672  l1rec->tg_sol[ipb] = exp(-A / l1rec->csolz[ip]);
673  l1rec->tg_sen[ipb] = exp(-A / l1rec->csenz[ip]);
674  }
675  }
676  } else {
677  if ((input->gas_opt & GAS_TRANS_TBL_BIT) != 0) {
678  if (firstRun) {
680  firstRun = FALSE;
681  }
682  /*
683  o2 transmittance is special
684  If the oxaband_opt is set to use the gas_transmittance tables, oblige
685  */
686  amf_solz=1.0/l1rec->csolz[ip];
687  amf_senz=1.0/l1rec->csenz[ip];
688  amf_total=amf_solz+amf_senz;
689 
690  if (amf) {
691  index_amf_senz = get_index_lowerbound(amf_mixed, num_airmass, amf_senz);
692  index_amf_solz = get_index_lowerbound(amf_mixed, num_airmass, amf_solz);
693  index_amf_total = get_index_lowerbound(amf_mixed, num_airmass, amf_total);
694 
695  ratio_senz = (amf_senz - amf_mixed[index_amf_senz]) /
696  (amf_mixed[index_amf_senz + 1] - amf_mixed[index_amf_senz]);
697  ratio_solz = (amf_solz - amf_mixed[index_amf_solz]) /
698  (amf_mixed[index_amf_solz + 1] - amf_mixed[index_amf_solz]);
699  ratio_total = (amf_total - amf_mixed[index_amf_total]) /
700  (amf_mixed[index_amf_total + 1] - amf_mixed[index_amf_total]);
701  }
702 
703  if (input->oxaband_opt == 2) {
704  o2_transmittance(l1rec, ip);
705  }
706  }
707 
708  if ((input->gas_opt & O3_BIT) != 0) {
710  }
711 
712  if ((input->gas_opt & CO2_BIT) != 0) {
714  }
715 
716  if ((input->gas_opt & NO2_BIT) != 0) {
718  }
719 
720  if ((input->gas_opt & H2O_BIT) != 0) {
722  }
723 
724  if ((input->gas_opt & CO_BIT) != 0) {
725  co_transmittance(l1rec, ip);
726  }
727  if ((input->gas_opt & CH4_BIT) != 0) {
729  }
730  if ((input->gas_opt & N2O_BIT) != 0) {
732  }
733  if (amf) {
734  for (ib = 0; ib < nwave; ib++) {
735  ipb = ip * nwave + ib;
736  l1rec->tg_sen[ipb] = l1rec->tg[ipb] / l1rec->tg_sol[ipb];
737  }
738  }
739  }
740 }
741 
743  filehandle *l1file = l1rec->l1file;
744  uncertainty_t *uncertainty=l1rec->uncertainty;
745  int32_t nwave = l1file->nbands;
746  int npix = l1rec->npix;
747  int32_t ipb, ip, ib;
748  float tg_oz, tg_co2, tg_no2, dt_oz, dt_co2, dt_no2;
749 
750  float mu0, mu, tau_oz;
751 
752  float a_285, a_225;
753  float tau_to200, dtau_to200;
754  float no2_tr200, dno2_tr200;
755 
756  if(!t_co2)
758 
759 
760  for (ip = 0; ip < npix; ip++) {
761  mu0 = l1rec->csolz[ip];
762  mu = l1rec->csenz[ip];
763 
764  if (l1rec->no2_tropo[ip] > 0.0) {
765  no2_tr200 = l1rec->no2_frac[ip] * l1rec->no2_tropo[ip];
766  dno2_tr200 = l1rec->no2_frac[ip] * uncertainty->dno2_tropo[ip];
767  } else {
768  no2_tr200 = 0.0;
769  dno2_tr200 = 0.0;
770  }
771 
772  for (ib = 0; ib < nwave; ib++) {
773  ipb = ip * nwave + ib;
774 
775  tau_oz = l1rec->oz[ip] * l1file->k_oz[ib];
776 
777  /* calculate error in tg_sol */
778  tg_oz = exp(-tau_oz / mu0);
779  dt_oz = tg_oz / mu0 * l1file->k_oz[ib] * uncertainty->doz[ip];
780 
781  tg_co2 = pow(t_co2[ib], 1.0 / mu0);
782  dt_co2 = 0.0;
783 
784  tg_no2 = 1.0;
785  dt_no2 = 0.0;
786 
787  if (l1file->k_no2[ib] > 0) {
788  a_285 = l1file->k_no2[ib] * (1.0 - 0.003 * (285.0 - 294.0));
789  a_225 = l1file->k_no2[ib] * (1.0 - 0.003 * (225.0 - 294.0));
790 
791  tau_to200 = a_285 * no2_tr200 + a_225 * l1rec->no2_strat[ip];
792  dtau_to200 = sqrt(pow(a_285*dno2_tr200, 2) + pow(a_225 * uncertainty->dno2_strat[ip], 2));
793 
794  tg_no2 = exp(-tau_to200 / mu0);
795  dt_no2 = tg_no2 / mu0*dtau_to200;
796  }
797  uncertainty->dtg_sol[ipb] = sqrt(pow(tg_no2 * tg_co2*dt_oz, 2) + pow(tg_oz * tg_co2*dt_no2, 2) + pow(tg_no2 * tg_oz*dt_co2, 2));
798 
799 
800  /* calculate error in tg_sen */
801  tg_oz = exp(-tau_oz / mu);
802  dt_oz = tg_oz / mu * l1file->k_oz[ib] * uncertainty->doz[ip];
803 
804  tg_co2 = pow(t_co2[ib], 1.0 / mu);
805  dt_co2 = 0.0;
806 
807  tg_no2 = 1.0;
808  dt_no2 = 0.0;
809  if (l1file->k_no2[ib] > 0) {
810  a_285 = l1file->k_no2[ib] * (1.0 - 0.003 * (285.0 - 294.0));
811  a_225 = l1file->k_no2[ib] * (1.0 - 0.003 * (225.0 - 294.0));
812 
813  tau_to200 = a_285 * no2_tr200 + a_225 * l1rec->no2_strat[ip];
814 
815  dtau_to200 = sqrt(pow(a_285*dno2_tr200, 2) + pow(a_225 * uncertainty->dno2_strat[ip], 2));
816 
817  tg_no2 = exp(-tau_to200 / mu);
818  dt_no2 = tg_no2 / mu*dtau_to200;
819  }
820  uncertainty->dtg_sen[ipb] = sqrt(pow(tg_no2 * tg_co2*dt_oz, 2) + pow(tg_oz * tg_co2*dt_no2, 2) + pow(tg_no2 * tg_oz*dt_co2, 2));
821  }
822  }
823 }
824 
825 float get_wv_band_ratio(l1str *l1rec,int32_t ip,float window1, float absorp_band,float window2){
826 
827  float wv;
828  int i;
829  static int firstcall=1;
830 
831  filehandle *l1file = l1rec->l1file;
832  int32_t nwave = l1file->nbands;
833  int32_t ipb = ip*nwave,index;
834  float * wave=l1file->fwave;
835  float *Lt=&l1rec->Lt[ipb];
836  // float u=l1rec->csenz[ip];
837  float u0=l1rec->csolz[ip];
838  float *F0=l1rec->Fo;
839  float rhot_interp,trans_wv_true,ratio_wv_total;
840  static float *rhot, *tran_interp;
841  int index_amf_wv_total;
842 
843 
844  int band1,band2,band_absorp;
845 
846  if(firstcall){
847  firstcall=0;
848  rhot=(float *)malloc(nwave*sizeof(float));
849  tran_interp=(float *)malloc(num_water_vapors*sizeof(float));
850  }
851 
852  band1=windex(window1,wave,nwave);
853  band2=windex(window2,wave,nwave);
854  band_absorp=windex(absorp_band,wave,nwave);
855 
856  for(i=band1;i<=band2;i++)
857  rhot[i]=PI*Lt[i]/F0[i]/u0;
858 
859  rhot_interp=rhot[band1]+(absorp_band-window1)*(rhot[band2]-rhot[band1])/(window2-window1);
860 
861  trans_wv_true=rhot[band_absorp]/rhot_interp;
862 
863  if (amf) {
864  index_amf_wv_total = get_index_lowerbound(amf_wv, num_airmass, amf_total);
865  ratio_wv_total = (amf_total - amf_wv[index_amf_wv_total]) /(amf_wv[index_amf_wv_total + 1] - amf_wv[index_amf_wv_total]);
866  index = (model * num_wavelengths + band_absorp) * num_airmass * num_water_vapors +index_amf_wv_total * num_water_vapors;
867  for (i = 0; i < num_water_vapors; i++)
868  tran_interp[i] = wvtbl[index + i] * ratio_wv_total +
869  wvtbl[index + num_water_vapors + i] * (1 - ratio_wv_total);
870 
871  for (i = 0; i < num_water_vapors; i++) {
872  if (trans_wv_true >= tran_interp[i])
873  break;
874  }
875  if (i == 0)
876  i = 1;
877  if (i == num_water_vapors)
878  i = num_water_vapors - 1;
879 
880  wv = cwv_all[i] + (trans_wv_true - tran_interp[i]) * (cwv_all[i] - cwv_all[i - 1]) /
881  (tran_interp[i] - tran_interp[i - 1]);
882  }else{
883  wv=l1rec->wv[ip];
884  }
885 
886  return (wv);
887 }
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
#define MAX(A, B)
Definition: swl0_utils.h:25
#define MIN(x, y)
Definition: rice.h:169
Utility functions for allocating and freeing three-dimensional arrays of various types.
#define H2O_BIT
Definition: l12_parms.h:60
size_t num_wavelengths
Definition: gas_trans.c:22
#define N2O_BIT
Definition: l12_parms.h:65
size_t num_airmass
Definition: gas_trans.c:22
#define FALSE
Definition: rice.h:164
#define NULL
Definition: decode_rs.h:63
map< string, float > F0
Definition: DDSensor.cpp:39
read l1rec
#define CO2_BIT
Definition: l12_parms.h:58
float mu
#define TRUE
Definition: rice.h:165
#define O3_BIT
Definition: l12_parms.h:57
instr * input
#define M_PI
Definition: dtranbrdf.cpp:19
#define PI
Definition: l3_get_org.c:6
const float A
Definition: calc_par.c:100
void load_gas_tables(l1str *l1rec)
Definition: gas_trans.c:38
int get_atrem_cor(l1str *l1rec, int32_t ip, float *rhot, float *tg_tot, float *tg_sol, float *tg_sen)
void no2_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:471
#define ATREM_BIT
Definition: l12_parms.h:61
int32_t model
Definition: gas_trans.c:20
void h2o_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:510
void gaseous_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:633
#define CH4_BIT
Definition: l12_parms.h:64
void ch4_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:375
#define CO_BIT
Definition: l12_parms.h:63
l1_input_t * l1_input
Definition: l1_options.c:9
#define NO2_BIT
Definition: l12_parms.h:59
void o2_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:408
float ja
Definition: atrem_cor.h:114
void co_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:340
void co2_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:304
#define ATREM_O2
Definition: atrem_corl1.h:27
size_t num_models
Definition: gas_trans.c:22
int32_t get_index_lowerbound(float *table_val, int32_t num_val, float val)
Definition: gas_trans.c:275
size_t num_water_vapors
Definition: gas_trans.c:22
float get_wv_band_ratio(l1str *l1rec, int32_t ip, float window1, float absorp_band, float window2)
Definition: gas_trans.c:825
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
#define GAS_TRANS_TBL_BIT
Definition: l12_parms.h:62
float mu0
void ozone_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:288
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
void gas_trans_uncertainty(l1str *l1rec)
Definition: gas_trans.c:742
int i
Definition: decode_rs.h:71
int npix
Definition: get_cmp.c:28
void n2o_transmittance(l1str *l1rec, int32_t ip)
Definition: gas_trans.c:438