ocssw V2020
get_atrem_corl1.c
Go to the documentation of this file.
1 /*
2  * get_atrem_cor2.c
3  *
4  * Created on: Feb 19, 2015
5  * Author: Rick Healy (richard.healy@nasa.gov)
6  *
7  * * Notes about water vapor VMRS and related quantities: *
8  * *
9  * VAPVRT(60) - a table containing 60 column vapor values (in unit of cm) *
10  * *
11  * VAP_SLANT(I) = VAPVRT(I) * 2.0, VAP_SLANT is a new table for containing *
12  * two-way total vapor amounts. Here the number "2" can be *
13  * changed to other numbers, e.g., 2.5, without major *
14  * effects on retrieved water vapor values. *
15  * *
16  * G_VAP(I = 1,..., NL) = true vapor geometric factor for each layer in *
17  * the model atmosphere (after adjusting for the elevated *
18  * surface. *
19  * *
20  * VMRM(I) = VMRM(I)*G_VAP(I). The VMRS are multiplied by the geometrical *
21  * factor. We can calculate the vapor transmittance on the *
22  * Sun-surface-sensor path by assuming a vertical path in *
23  * the model atmosphere with geometric-factor-adjusted VMRS. *
24  * *
25  * CLMVAP = vertical column amount from ground to space in model atmosphere*
26  * CLMVAPP = vertical column amount from ground to aircraft or satellite *
27  * sensor in model atmosphere *
28  * Q = 2.152E25 = # of molecules above the surface at one atmosphere *
29  * (in unit of molecules/cm**2) *
30  * *
31  * VAP_SLANT_MDL= CLMVAP/COS(SOLZNI) + CLMVAPP/COS(OBSZNI) = total amount *
32  * of water vapor in the model atmosphere in the L-shaped *
33  * Sun-surface-plane ray path. *
34  * *
35  * G_VAP_EQUIV = VAP_SLANT_MDL / CLMVAP = the "equivalent" geometrical *
36  * factor corresponding to the total slant vapor amount *
37  * VAP_SLANT_MDL and the column vapor amount CLMVAP. *
38  * *
39  * SSH2O(I) (I = 1, ..., 60) - a pure scaling factor relative to the total *
40  * slant vapor amount of VAP_SLANT_MDL, and *
41  * SSH2O(I) = VAP_SLANT(I) / VAP_SLANT_MDL *
42  * *
43  * SH2O = one value of SSH2O(I). SH2O is used during generation of the *
44  * look-up table. *
45  * *
46  * VAPTT = VAP_SLANT_MDL*SH2O, is the absolute total vapor amount on the *
47  * L-shaped path corresponding to a spectrum stored in the *
48  * look-up table. *
49  * *
50  * CLMWVP = 0.5*(VAPTTA+VAPTTB)/G_VAP_EQUIV, is the retrieved column water *
51  * vapor amount from imaging spectrometer data. *
52  ********************************************************************************
53  *
54  *
55  * 01/15/2016 - Split the solar and sensor paths directly inside atrem.
56  * specify with atrem_splitpaths. Only works with atrem_full on.
57  *
58  * 10/20/2015 - r.healy - added ATREM options to command line
59  * atrem_opt - select gases for transmittance calculation
60  * atrem_geom - turn on/off geometry recalculation
61  * 0 = geometry calculated on error angle limit
62  * 1 = recalculate every pixel
63  * atrem_full - turn on/off full calculation (k-dist otherwise)
64  * atrem_model - select atmospheric model (1-6)
65  * - 0= determine from latitude and day of year
66  *
67  * !* References: *
68 !* Gao, B.-C., K. Heidebrecht, and A. F. H. Goetz, Derivation of scaled *
69 !* surface reflectances from AVIRIS data, Remote Sens. Env., 44, *
70 !* 165-178, 1993. *
71 !* Gao, B.-C., and C. O. Davis, Development of an operational algorithm *
72 !* for removing atmospheric effects from HYDICE and HSI data, *
73 !* in SPIE'96 Conference Proceedings, Vol. 2819, 45-55, 1996. *
74 !* Gao, B.-C., and A. F. H. Goetz, Column atmospheric water vapor and *
75 !* vegetation liquid water retrievals from airborne imaging *
76 !* spectrometer data, J. Geophys. Res., 95, 3549-3564, 1990. *
77 !* Goetz, A. F. H., and M. Herring, The high resolution imaging *
78 !* spectrometer (HIRIS) for Eos, IEEE Trans. Geosci. Remote Sens., 27,*
79 !* 136-144, 1989. *
80 !* Goetz, A. F. H., G. Vane, J. Solomon, and B. N. Rock, Imaging *
81 !* spectrometry for Earth remote sensing, Science, 228, 1147-1153,1985*
82 !* Green, R. O., and B.-C. Gao, A Proposed Update to the Solar Irradiance*
83 !* Spectrum used in LOWTRAN and MODTRAN, in Summaries of the Fourth *
84 !* Annual JPL Airborne Geoscience Workshop, October 25-29, (Editor, *
85 !* R. O. Green), JPL Publ. 93-26, Vol. 1, pp. 81-84, Jet Propul. Lab, *
86 !* Pasadena, Calif., 1993. *
87 !* Kneizys, F. X., E. P. Shettle, L. W. Abreu, J. H. Chetwynd, G. P. *
88 !* Anderson, W. O. Gallery, J. E. A. Selby, and S. A. Clough, Users *
89 !* guide to LOWTRAN7, AFGL-TR-8-0177, Air Force Geophys. Lab., *
90 !* Bedford, Mass., 1988. *
91 !* Iqbal, M., An Introduction To Solar Radiation, Academic, San Diego, *
92 !* Calif., 1983. *
93 !* Malkmus, W., Random Lorentz band model with exponential-tailed S line *
94 !* intensity distribution function, J. Opt. Soc. Am., 57, 323-329,1967*
95 !* Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, *
96 !* Numerical Recipes-The ART of Scientific Computing, Cambridge *
97 !* University Press, 1986. *
98 !* Rothman, L. S., et al., The HITRAN 2008 molecular spectroscopic *
99 !* database, JQSRT, 110, 533-572, 2009. *
100 !* Solomon, S., R. W. Portmann, R. W. Sanders, J. S. Daniel, W. Madsen, *
101 !* B. Bartram, and E. G. Dutton, On the role of nitrogen dioxide in *
102 !* the absorption of solar radiation, J. Geophys. Res., 104, *
103 !* 12047-12058, 1999. *
104 !* Tanre, D., C. Deroo, P. Duhaut, M. Herman, J. J. Morcrette, J. Perbos,*
105 !* and P. Y. Deschamps, Description of a computer code to simulate *
106 !* the satellite signal in the solar spectrum: the 5S code, Int. *
107 !* J. Remote Sens., 11, 659-668, 1990. *
108 !* Tanre, D., C. Deroo, P. Duhaut, M. Herman, J. J. Morcrette, J. Perbos,*
109 !* and P. Y. Deschamps, Simulation of the satellite signal in the *
110 !* solar spectrum (5S), Users' Guide (U. S. T. De Lille, 59655 *
111 !* Villeneu d'Ascq, France: Laboratoire d'Optique Atmospherique), *
112 !* 1986. *
113 !* Thuillier, G., et al., Solar irradiance reference spectra for two *
114 !* solar active levels, Adv. Space Res., 34, 256-261, 2004. *
115 !* Vane, G., R. O. Green, T. G. Chrien, H. T. Enmark, E. G. Hansen, and *
116 !* W. M. Porter, The Airborne Visible/Infrared Imaging Spectrometer, *
117 !* Remote Sens. Env., 44, 127-143, 1993. *
118 !* Vane, G. (Ed), Airborne visible/infrared imaging spectrometer *
119 !* (AVIRIS), JPL Publ. 87-38, Jet Propul. Lab, Pasadena, Calif., 1987.*
120 !* Vermote, E., D. Tanre, J. L. Deuze, M. Herman, and J. J. Morcrette, *
121 !* Second simulation of the satellite signal in the solar spectrum *
122 !* (6S), 6S User's Guide Version 1, NASA-GSFC, Greenbelt, Maryland, *
123 !* 134 pages, 1994. *
124 !* *
125 !********************************************************************************
126  *
127  */
128
129 /* compile this code with:
130  gcc -O1 -Wpadded -Wpacked -malign-double -mpreferred-stack-boundary=8 -o get_atrem_cor3 \
131  get_atrem_cor3.c rdatreminfo.c atrem_app_refl_plus_gas_removal_for_l2gen3.o cubeio.o \
132  tpvmr_init.o solar_irr_PC.o bndprms.o -lgfortran \
133  -L/disk01/home/rhealy/ocssw/build/cbuild/src/libgenutils -lgenutils \
134  -L/disk01/home/rhealy/ocssw/build/lib3/src/netcdf/netcdf-fortran-4.2/fortran/.libs \
135  -lnetcdff -L/disk01/home/rhealy/ocssw/build/lib3/src/netcdf/netcdf-4.3.1.1/liblib/.libs -lnetcdf \
136  -L/disk01/home/rhealy/ocssw/build/lib3/src/hdf5/hdf5-1.8.10-patch1/src/.libs -lhdf5 \
137  -L/disk01/home/rhealy/ocssw/build/lib3/src/hdf5/hdf5-1.8.10-patch1/hl/src/.libs -lhdf5_hl \
138  -L/disk01/home/rhealy/ocssw/build/lib3/src/netcdf/netcdf-4.3.1.1/libsrc4/.libs -lnetcdf4 -lz \
139  -L/disk01/home/rhealy/ocssw/build/lib3/src/hdf5/hdf5-1.8.10-patch1/src/.libs -lhdf5
140
141  * compile fortran routines:
142  gfortran -malign-double -c atrem_app_refl_plus_gas_removal_for_l2gen2.f90 cubeio.f90 \
143  bndprms.f solar_irr_PC.f tpvmr_init.f
144
145 ---- Previous version without netcdf -----
146
147  gcc -Wpadded -Wpacked -malign-double -mpreferred-stack-boundary=8 -o get_atrem_cor \
148  get_atrem_cor.c atrem_app_refl_plus_gas_removal_for_l2gen.o cubeio.o \
149  tpvmr_init.o solar_irr_PC.o bndprms.o -lgfortran -L/disk01/home/rhealy/ocssw/build/cbuild/src/libgenutils -lgenutils
150  *
151  * compile fortran routines:
152  gfortran -malign-double -I/disk01/home/rhealy/ocssw/build/lib3/src/netcdf/netcdf-fortran-4.2/fortran -c atrem_app_refl_plus_gas_removal_for_l2gen3.f90 cubeio.f90 \
153  bndprms.f solar_irr_PC.f tpvmr_init.f
154
155  Example:
156  get_atrem_cor2 < input/test_input.txt > cor2_test90b.out
157  */
158 #include "atrem_corl1.h"
159 #include <sensorDefs.h>
160 #include <timeutils.h>
161 #include <math.h>
162 #define BUFSZ 100
163 #define MINDEGCHANGE 55 // Minimum change (degrees squared) in sum of zenith and azimuth angle squared
164 // for both solar and sensor angles before recalculating tran_table (transmittance table)
165
166 int get_atrem_cor(int32_t sensorID, l1str *l1rec, int32_t ip, float *rhot, float *tg_tot, float *tg_sol, float *tg_sen) {
167
168  static paramstr P;
169  int i, j, nb, modnum = input->atrem_model;
170  static int firstCall = 1;
171  int32_t nbands = l1rec->l1file->nbands;
172  static double prev_ddeg = -999, prev_max_senz, prev_min_senz;
173  static float **angle_limit, *ang_senz, *ang_solz;
174  static int n_senz, n_solz, so_ang, se_ang;
175  static int prev_modnum;
176  float limitang, *anglelimit;
177
178  //Initialize input paramters
179
180
181  if (firstCall == 1) { //|| prev_dist > MINDEGCHANGE) {
182  //INIT
184  prev_modnum = P.model;
185
186  if (P.dogeom == 0) {
187  printf("Reading geometry limiting angles for Atrem\n");
188  if (get_angle_limits(&anglelimit, &ang_senz, &ang_solz, &n_senz, &n_solz)) {
189  printf("-E- %s line %d : Error reading angle_limit file.\n",
190  __FILE__, __LINE__);
191  exit(FATAL_ERROR);
192  }
193
194  if ((angle_limit = (float **) calloc(n_senz, sizeof (float *))) == NULL) {
195  printf("-E- : Error allocating memory to tindx\n");
196  exit(FATAL_ERROR);
197  }
198
199  for (j = 0; j < n_senz; j++) {
200  if ((angle_limit[j] = (float *) calloc(n_solz, sizeof (float))) == NULL) {
201  printf("-E- : Error allocating memory to tindx\n");
202  exit(FATAL_ERROR);
203  }
204  for (i = 0; i < n_solz; i++) {
205  angle_limit[j][i] = *(anglelimit + j * (n_solz) + i);
206  //printf("RJH:2: angle_limit: %f %f %f\n",ang_solz[i],ang_senz[j],angle_limit[j][i]);
207
208  }
209  }
210
211  prev_max_senz = l1rec->senz[ip];
212  prev_min_senz = l1rec->senz[ip];
213  }
214  }
215
216  if (P.dogeom == 0) {
217  if (l1rec->senz[ip] > prev_max_senz) prev_max_senz = l1rec->senz[ip];
218  if (l1rec->senz[ip] < prev_min_senz) prev_min_senz = l1rec->senz[ip];
219
220  prev_ddeg = fabs(prev_max_senz - prev_min_senz);
221
222  limitang = get_current_angle_limit(l1rec->senz[ip], l1rec->solz[ip], &se_ang, &so_ang, angle_limit, ang_senz, ang_solz, n_senz, n_solz);
223  //printf("RJH: prev_max_senz=%f prev_min_senz=%f limitang[%d]=%f prev_ddeg=%f senz=%f solz=%f \n",prev_max_senz,prev_min_senz,ip,limitang,prev_ddeg,l1rec->senz[ip],l1rec->solz[ip]);
224  }
225  if (firstCall == 1 || prev_ddeg >= limitang || P.dogeom != 0) {
226
227  //printf("Calculating Transmittance table for Atrem correction, ip=%d\n",ip);
228  //start_time = now();
229  //printf("\nBegin get_atrem_cor processing at %s\n\n", ydhmsf(start_time,'L'));
230
231
232  geometry_l2gen_.splitpaths = 0;
233  geometry_l2gen_.senzn_l2 = l1rec->senz[ip];
234  geometry_l2gen_.senaz_l2 = l1rec->sena[ip];
235  geometry_l2gen_.solzn_l2 = l1rec->solz[ip];
236  getinput3_.vrto3 = l1rec->oz[ip];
237
238  geometry_();
239  //printf("Processed geometry after %6.0f seconds\n",now()-start_time);
240
241  init_speccal_();
242  //printf("Processed init_speccal after %6.0f seconds\n",now()-start_time);
243
244  tran_table_();
245  //printf("Processed tran_table after %6.0f seconds\n",now()-start_time);
246
247  P.nh2o = init_speccal3_.nh2o;
248  P.start_ndx[0] = init_speccal6_.ist1 - 1; // Fortran array index starts at 1, C at 0
249  P.start_ndx[1] = init_speccal6_.ist2 - 1;
250  P.end_ndx[0] = init_speccal6_.ied1 - 1;
251  P.end_ndx[1] = init_speccal6_.ied2 - 1;
252  P.start_p94 = init_speccal6_.istp94 - 1;
253  P.end_p94 = init_speccal6_.iedp94 - 1;
254  P.start_ndx[2] = init_speccal7_.ist3 - 1;
255  P.start_ndx[3] = init_speccal7_.ist4 - 1;
256  P.end_ndx[2] = init_speccal7_.ied3 - 1;
257  P.end_ndx[3] = init_speccal7_.ied4 - 1;
258  P.start_1p14 = init_speccal7_.ist1p14 - 1;
259  P.end_1p14 = init_speccal7_.ied1p14 - 1;
260  P.wt1 = init_speccal8_.wt1;
261  P.wt2 = init_speccal8_.wt2;
262  P.wt3 = init_speccal8_.wt3;
263  P.wt4 = init_speccal8_.wt4;
264  P.start2 = init_speccal10_.istrt2 - 1;
265  P.end2 = init_speccal10_.iend2 - 1;
266  P.ncv2 = init_speccal10_.ncv2;
267  P.finst2 = init_speccal10_.finst2;
268  P.natot = init_speccal11_.natot;
269  P.nbtot = init_speccal11_.nbtot;
270  P.nctot = init_speccal11_.nctot;
271  P.ndtot = init_speccal11_.ndtot;
272  P.g_vap_equiv = geometry3_.g_vap_equiv;
273  P.r0p94 = tran_table1_.r0p94;
274  P.r1p14 = tran_table1_.r1p14;
275  P.vaptot = tran_table1_.vaptot;
276  P.trntbl = tran_table1_.trntbl;
277
278
279  // printf("nb: %d %d %d %d\n",P.nb1,P.nb2,P.nb3,P.nb4);
280  // printf("nb: %d %d\n",P.nbp94,P.nb1p14);
281  // printf("nh2o: %d\n",P.nh2o);
282  // printf("nobs: %d\n",P.nobs);
283  // printf("startndx: %d %d %d %d\n",P.start_ndx[0],P.start_ndx[1],P.start_ndx[2],P.start_ndx[3]);
284  // printf("endndx: %d %d %d %d\n",P.end_ndx[0],P.end_ndx[1],P.end_ndx[2],P.end_ndx[3]);
285  // printf("start_p94: %d\n",P.start_p94);
286  // printf("end_p94: %d\n",P.end_p94);
287  // printf("start_1p14: %d\n",P.start_1p14);
288  // printf("end_1p14: %d\n",P.end_1p14);
289  // printf("%d\n",P.start2);
290  // printf("%d\n",P.end2);
291  //
292  // printf("ncv2: %d\n",P.ncv2);
293  // printf("ntot: %d %d %d %d\n",P.natot,P.nbtot,P.nctot,P.ndtot);
294  // printf("wt: %f %f %f %f\n",P.wt1,P.wt2,P.wt3,P.wt4);
295  // printf("delta: %f %f\n",P.delta,P.delta2);
296  // printf("g_vap_equiv: %f \n",P.g_vap_equiv);
297  // printf("ip=%d\n",ip);
298  // printf("prev_dist=%f\n",prev_ddeg);
299  firstCall = 0;
300
301  prev_max_senz = l1rec->senz[ip];
302  prev_min_senz = l1rec->senz[ip];
303  }
304
305  // if (l1rec->iscan != prevscan) {
306  // prevscan = l1rec->iscan;
307  // }
308  if (modnum == 0) {
309  int16_t year, day;
310  double secs;
311  unix2yds(l1rec->scantime, &year, &day, &secs);
312  P.model = getModelNum(*l1rec->lat, day);
313  if (P.model != prev_modnum) {
314  nb = init_tpvmr(P.model);
315  if (nb <= 0) {
316  printf("-E- %s line %d : Atmospheric data could not be initialized.\n",
317  __FILE__, __LINE__);
318  exit(FATAL_ERROR);
319  }
321  prev_modnum = P.model;
322  }
323  }
324  l1rec->wv[ip] = get_atrem(tg_tot, rhot, &P); //returns water vapor(cm)
325
326  // Calculate the separate path transmittances (solar and sensor) if selected
327  if (input->atrem_splitpaths) {
328  geometry_l2gen_.splitpaths = 1;
329  geometry_l2gen_.water_vapor = l1rec->wv[ip]; //Misunderstanding. This isn't used.
330  geometry_l2gen_.ja = P.ja;
331  geometry_l2gen_.jb = P.jb;
332  geometry_l2gen_.f1a = P.f1a;
333  geometry_l2gen_.f1b = P.f1b;
334  geometry_l2gen_.f2a = P.f2a;
335  geometry_l2gen_.f2b = P.f2b;
336
337  // Call geometry with the indices from the tran_table
338  geometry_();
339
340  // Don't need to call init_speccal since other trace gas transmittances don't depend on water vapor
341
342  // recalculate transmittances based on separate solar/sensor paths using the tran_table indices from the first call
343  tran_table_();
344
345  for (i = 0; i < P.nobs; i++) {
346  tg_sol[i] = tran_table_l2gen_.tg_sol[i];
347  tg_sen[i] = tran_table_l2gen_.tg_sen[i];
348  }
349
350  }
351  //printf("Water vapor[%d]=%f\n",ip,l1rec->wv[ip]);
352  //printf("Processed get_atrem after %6.0f seconds\n",now()-start_time);
353
354
355  return (0);
356 }
357
358 float get_atrem(float *tg_tot, float *rhot, paramstr *P) {
359
360  double const1 = 0, const2 = 0, const3 = 0, const4 = 0, const5 = 0, const6 = 0;
361  double ratio_94c, ratio_94co, ratio_114c, ratio_114co;
362  int32_t i, ja, jb;
363  float clmwvp;
364  /*
365  Arrays related to look-up table:
366  VAPTOT: TOTAL SUN-SURFACE-SENSOR PATH WATER VAPOR IN UNITS OF CM
367  R0P94 : Tc(0.94 um)/(WT1*Tc(0.86)+WT2*Tc(1.02))
368  R1P14 : Tc(1.14 um)/(WT3*Tc(1.05)+WT4*Tc(1.23))
369
370  Calculating 3-channel ratios from an observed spectrum, using a
371  look-up table procedure to derive the column amount of water vapor
372  and to find the associated simulated transmittance spectrum.
373  */
374  for (i = P->start_ndx[0]; i <= P->end_ndx[0]; i++) {
375  const1 += rhot[i];
376  }
377  const1 /= P->nb1;
378
379  for (i = P->start_ndx[1]; i <= P->end_ndx[1]; i++) {
380  const2 += rhot[i];
381  }
382  const2 /= P->nb2;
383
384  // printf("const2=%f nb2=%d istr2=%d ind2=%d\n",const2,P->nb2,P->start_ndx[1],P->end_ndx[1]);
385
386  for (i = P->start_p94; i <= P->end_p94; i++) {
387  const3 += rhot[i];
388  }
389  const3 /= P->nbp94;
390
391  ratio_94co = const3 / ((P->wt1 * const1) + (P->wt2 * const2));
392  ratio_94c = ratio_94co;
393
394  if (ratio_94co > 1.0) {
395  const1 = 0.0;
396
397  for (i = P->start_ndx[0]; i <= P->end_ndx[0]; i++) {
398  const1 += (1.0 / rhot[i]);
399  }
400  const1 /= P->nb1;
401
402  const2 = 0.0;
403  for (i = P->start_ndx[1]; i <= P->end_ndx[1]; i++) {
404  const2 += (1.0 / rhot[i]);
405  }
406  const2 /= P->nb2;
407  const3 = 0.0;
408  for (i = P->start_p94; i <= P->end_p94; i++) {
409  const3 += (1.0 / rhot[i]);
410  }
411  const3 /= P->nbp94;
412
413  ratio_94c = const3 / ((P->wt1 * const1) + (P->wt2 * const2));
414  }
415
416  debug_atrem.rp94 = ratio_94c;
417
418  const4 = 0.0;
419  for (i = P->start_ndx[2]; i <= P->end_ndx[2]; i++) {
420  const4 += rhot[i];
421  }
422  const4 /= P->nb3;
423
424  const5 = 0.0;
425  for (i = P->start_ndx[3]; i <= P->end_ndx[3]; i++) {
426  const5 += rhot[i];
427  }
428  const5 /= P->nb4;
429
430  const6 = 0.0;
431  for (i = P->start_1p14; i <= P->end_1p14; i++) {
432  const6 += rhot[i];
433  }
434  const6 /= P->nb1p14;
435
436  /* DEBUG
437  *
438  */
439  debug_atrem.cst1 = const1;
440  debug_atrem.cst2 = const2;
441  debug_atrem.cst3 = const3;
442  debug_atrem.cst4 = const4;
443  debug_atrem.cst5 = const5;
444  debug_atrem.cst6 = const6;
445
446  ratio_114co = const6 / ((P->wt3 * const4) + (P->wt4 * const5));
447  ratio_114c = ratio_114co;
448
449  if (ratio_114co > 1.0) {
450
451  const4 = 0.0;
452  for (i = P->start_ndx[2]; i <= P->end_ndx[2]; i++) {
453  const4 += (1.0 / rhot[i]);
454  }
455  const4 /= P->nb3;
456  for (i = P->start_ndx[3]; i <= P->end_ndx[3]; i++) {
457  const5 += (1.0 / rhot[i]);
458  }
459  const5 /= P->nb4;
460  const6 = 0.0;
461  for (i = P->start_1p14; i <= P->end_1p14; i++) {
462  const6 += (1.0 / rhot[i]);
463  }
464  const6 /= P->nb1p14;
465  ratio_114c = const6 / ((P->wt3 * const4) + (P->wt4 * const5));
466  }
467
468  debug_atrem.r1p14 = ratio_114c;
469
470  double delta, deltab, fja, fjap1, fjb, fjbp1, vaptta, vapttb;
471  double speca[NBANDS], specb[NBANDS], specav, spec450;
472  ja = P->nh2o / 2;
473  ja = hunt(P->r0p94, P->nh2o, ratio_94c, ja);
474  // printf("xx[%d]=%f xx[%d]=%f JA=%d\n",ja-1,P->r0p94[ja-1],ja,P->r0p94[ja],ja);
475  if (ja >= 0 && ja < NH2OMAX) {
476  delta = P->r0p94[ja + 1] - P->r0p94[ja];
477  fja = (P->r0p94[ja + 1] - ratio_94c) / delta;
478  fjap1 = (ratio_94c - P->r0p94[ja]) / delta;
479  vaptta = fja * P->vaptot[ja] + fjap1 * P->vaptot[ja + 1];
480  if (ratio_94co > 1.0) vaptta = -vaptta;
481  } else {
482  if (ja < 0) vaptta = P->vaptot[ja + 1];
483  if (ja > P->nh2o) vaptta = P->vaptot[ja];
484  }
485
486  if (ratio_94co <= 1.0) {
487  for (i = 0; i < P->nobs; i++) {
488  if (ja >= 0 && ja < NH2OMAXM1) {
489  speca[i] = fja * P->trntbl[ja][i] + fjap1 * P->trntbl[ja + 1][i];
490  } else {
491  if (ja < 0) speca[i] = P->trntbl[ja + 1][i];
492  if (ja >= NH2OMAXM1) speca[i] = P->trntbl[ja][i];
493  }
494  }
495  }
496
497  if (ratio_94co > 1.0) {
498  for (i = 0; i < P->nobs; i++) {
499  if (ja >= 0 && ja < NH2OMAXM1) {
500  speca[i] = 1.0 / (fja * P->trntbl[ja][i] + fjap1 * P->trntbl[ja + 1][i]);
501  } else {
502  if (ja < 0) speca[i] = 1.0 / P->trntbl[ja + 1][i];
503  if (ja >= NH2OMAXM1) speca[i] = 1.0 / P->trntbl[ja][i];
504  }
505  }
506  }
507
508  jb = ja;
509
510  jb = hunt(&P->r1p14[0], P->nh2o, ratio_114c, jb);
511  //printf("RJH: 1p14 JB=%d\n",jb);
512
513  debug_atrem.jac = ja;
514  debug_atrem.jbc = jb;
515
516  if (jb >= 0 && jb < NH2OMAXM1) {
517  deltab = P->r1p14[jb + 1] - P->r1p14[jb];
518  fjb = (P->r1p14[jb + 1] - ratio_114c) / deltab;
519  fjbp1 = (ratio_114c - P->r1p14[jb]) / deltab;
520  vapttb = fjb * P->vaptot[jb] + fjbp1 * P->vaptot[jb + 1];
521  if (ratio_114co > 1.0) vapttb = -vapttb;
522  } else {
523  if (jb < 0) vapttb = P->vaptot[jb + 1];
524  if (jb <= NH2OMAXM1) vapttb = P->vaptot[jb];
525  }
526
527  if (ratio_114co <= 1.0) {
528  for (i = 0; i < P->nobs; i++) {
529  if (jb >= 0 && jb < NH2OMAXM1) {
530  specb[i] = fjb * P->trntbl[jb][i] + fjbp1 * P->trntbl[jb + 1][i];
531  } else {
532  if (jb < 0) specb[i] = P->trntbl[jb + 1][i];
533  if (jb >= NH2OMAXM1) specb[i] = P->trntbl[jb][i];
534  }
535  }
536  }
537  if (ratio_114co > 1.0) {
538  for (i = 0; i < P->nobs; i++) {
539  if (jb >= 0 && jb < NH2OMAXM1) {
540  specb[i] = 1.0 / (fjb * P->trntbl[jb][i] + fjbp1 * P->trntbl[jb + 1][i]);
541  } else {
542  if (jb < 0) specb[i] = 1.0 / P->trntbl[jb + 1][i];
543  if (jb >= NH2OMAXM1) specb[i] = 1.0 / P->trntbl[jb][i];
544  }
545  }
546  }
547
548  clmwvp = 0.5 * (vaptta + vapttb) / P->g_vap_equiv;
549  spec450 = 1.0 - 0.5 * (speca[P->idx450] + specb[P->idx450]); // should be 0.0
550
551  // Derivation of surface reflectances
552
553  for (i = 0; i < P->nobs; i++) {
554  specav = 0.5 * (speca[i] + specb[i]);
555  tg_tot[i] = specav + spec450;
556  //printf("RJH: Atrem: %d %f YY=%f %f %f\n",i,tg_tot[i],rhot[i],rhot[i]/specav,rhot[i]-rhot[i]/specav);
557  }
558
559  P->ja = ja;
560  P->jb = jb;
561  P->f1a = fja;
562  P->f2a = fjap1;
563  P->f1b = fjb;
564  P->f2b = fjbp1;
565
566
567  /*
568  // Smooth the derived surface reflectance spectra
569  // rhot = yy from atrem fortran code
570  // tg_tot = vc from atrem fortran code
571  int ii,j;
572  double truncv;
573
574  if (P->delta2 > P->delta) {
575
576  // * First, replace radiances <= 0 near the spectral overlapping parts of the
577  // * four AVIRIS spectrometers by radiances in the nearby AVIRIS' channels.
578
579  for (i=P->natot-3; i< P->natot+2; i++) {
580  if (rhot[i] <= 0.0) rhot[i] = rhot[i-1];
581  }
582  for (i=P->nbtot-3; i< P->nbtot+2; i++) {
583  if (rhot[i] <= 0.0) rhot[i] = rhot[i-1];
584  }
585  for (i=P->nctot-3; i< P->nctot+2; i++) {
586  if (rhot[i] <= 0.0) rhot[i] = rhot[i-1];
587  }
588  for (i=P->ndtot-3; i< P->ndtot+2; i++) {
589  if (rhot[i] <= 0.0) rhot[i] = rhot[i-1];
590  }
591  for (i=P->start2-1; i<P->end2; i++){
592  truncv = 0.0;
593  ii = i - P->ncv2 - 1;
594  for (j=ii; j<i+P->ncv2; j++) {
595  truncv += rhot[j]*P->finst2[j-ii+1];
596  }
597  rhot[i] = truncv;
598  }
599  }
600  */
601
602  return (clmwvp);
603
604 }
605
606 /*********************************************************************************
607  * *
608  * Name: HUNT *
609  * Purpose: finds the element in array XX that is closest to value X. Array AA *
610  * must be monotonic, either increasing or decreasing. *
611  * Parameters: XX - array to search *
612  * N - number of elements in the array *
613  * X - element to search for closest match *
614  * JLO - index of the closest matching element *
615  * Algorithm: this subroutine was copied from Numerical Recipes *
616  * Modified for C and cleaned up by *
617  * Richard Healy SAIC/NASA-GSFC 2/19/2015 *
618  * (Who didn't realize it was a Numerical Recipe's function until finished) *
619  * Globals used: none *
620  * Global output: none *
621  * Return Codes: none *
622  * Special Considerations: none *
623  * *
624  **********************************************************************************
625  *
626  */
627 int32_t hunt(float *xx, int32_t n, double x, int32_t jlo) {
628  int32_t inc, jhi, jm;
629  int ascnd;
630
631  ascnd = (xx[n - 1] >= xx[0]);
632
633  inc = 1;
634  if (jlo >= 0 && jlo < n) {
635
636  if ((x >= xx[jlo]) == ascnd) {
637  jhi = jlo + inc;
638  while (jhi < n && ((x >= xx[jhi]) == ascnd)) {
639  jlo = jhi;
640  inc += inc;
641  jhi = jlo + inc;
642  }
643  if (jhi >= n)
644  jhi = n;
645
646  } else {
647  jhi = jlo;
648  jlo = jhi - inc;
649  while (jlo >= 0 && ((x < xx[jlo]) == ascnd)) {
650  jhi = jlo;
651  inc += inc;
652  jlo = jhi - inc;
653  }
654  if (jlo < 0)
655  jlo = -1;
656
657  }
658
659  } else {
660  jlo = 0;
661  jhi = n + 1;
662  }
663  while (jhi - jlo != 1) {
664
665  jm = (jhi + jlo) / 2;
666
667  if ((x >= xx[jm]) == ascnd) {
668  jlo = jm;
669  } else {
670  jhi = jm;
671  }
672  }
673
674  if (x == xx[n - 1]) jlo = n - 2;
675  if (x == xx[0]) jlo = 0;
676
677  return (jlo);
678 }
679
680 int init_tpvmr(int model) {
681
682  int i, nb;
683  printf("RJH: Initializing ATM for model number = %d\n", model);
684  model--;
685  if (model < 0 || model > MODELMAX) {
686  printf("-E- %sline %d: Invalid ATM Model number\n Value must be between 1 and 7\n: get_atrem_cor3\n",
687  __FILE__, __LINE__);
688  exit(FATAL_ERROR);
689  }
690  nb = tpvmr_init1_.tpvmr[0][model];
691
692  for (i = 0; i < nb; i++) {
693  getinput3_.h[i] = tpvmr_init1_.tpvmr[1 + (4 * i)][model];
694  //Convert the atmospheric pressure from "mb" to "atm.":
695  getinput3_.p[i] = tpvmr_init1_.tpvmr[2 + (4 * i)][model] / 1013.;
696  getinput3_.t[i] = tpvmr_init1_.tpvmr[3 + (4 * i)][model];
697  //Convert the VMR from the ppm unit in the model to absolute unit
698  getinput3_.vmr[i] = tpvmr_init1_.tpvmr[4 + (4 * i)][model]*1.0E-06;
699  }
700
701  for (i = nb; i < MODELMAX; i++) {
702  getinput3_.h[i] = 1000.;
703  getinput3_.p[i] = 0.0;
704  getinput3_.t[i] = 300;
705  getinput3_.vmr[i] = 0.0;
706  }
707  getinput3_.nb = nb;
708  getinput3_.nl = nb - 1;
709
710
711
712  return nb;
713 }
714
715 int getModelNum(float lat, int32_t day) {
716  // Determine which atmospheric model to use
717  // 1 = tropical
718  // 2 = mid latitude summer
719  // 3 = mid latitude winter
720  // 4 = subarctic summer
721  // 5 = subarctic winter
722  // 6 = US standard 1962
723  int16 mon = (int) day / 31 + 1; // month of year (no need for perfection..at least according to the sea surface salinity reference algorithm)
724
725  if (fabs(lat) < 30) return (1);
726  else {
727  switch (mon) {
728  case 12:
729  case 1:
730  case 2:
731  if (lat < 60 && lat > 30) return (3);
732  if (lat < -30 && lat >-60) return (2);
733  if (lat > 60) return (5);
734  return (4);
735  break;
736  case 6:
737  case 7:
738  case 8:
739  if (lat < 60 && lat > 30) return (2);
740  if (lat < -30 && lat >-60) return (3);
741  if (lat > 60) return (4);
742  return (5);
743  break;
744  default:
745  return (6);
746
747  }
748
749
750  }
751
752
753 }
754
755 int init_atrem(int32_t sensorID, paramstr *P, l1str *l1rec, int32_t nbands) {
756
757  int32_t nb, atrem_opt = input->atrem_opt, atrem_full = input->atrem_full, atrem_geom = input->atrem_geom;
758  int32_t atrem_splitpaths = input->atrem_splitpaths, atrem_model = input->atrem_model, gas_opt = input->gas_opt;
759  float *fwhm, *xppp;
760  float nwave;
761  int i;
762  char *filedir;
763  char filename[FILENAME_MAX];
764  int *model, flag_gas;
765
766  model = (int *) calloc(1, sizeof (model));
767  fwhm = (float *) calloc(1, sizeof (fwhm));
768  xppp = (float *) calloc(1, sizeof (xppp));
769
770  if ((filedir = getenv("OCDATAROOT")) == NULL) {
771  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
772  return (-1);
773  }
774
775  strcpy(filename, filedir);
776  strcat(filename, "/");
778  strcat(filename, "/\0");
779
780  //INIT - get_input_ must be called before any common block structures are filled in.
781  // It will also fill some of the parameters with default values
782  get_input_();
783
784  strcpy(input_l2gen_.filename, filename);
785  input_l2gen_.dln = strlen(filename);
786  //printf("filedir=%s Filename=%s len=%d\n",filedir, filename,input_l2gen_.dln);
787
788  if (!(fwhm = l1rec->l1file->fwhm)) {
789  nwave = rdatreminfo(sensorID, 0, "fwhm", (void **) &fwhm);
790  }
791  for (i = 0; i < nbands; i++) {
792  getinput4_.wavobs[i] = l1rec->l1file->fwave[i] / 1000.; // Convert nm to microns
793  getinput4_.fwhm[i] = fwhm[i];
794  //printf("->RJH:fwave(%d)=%f fwhm(%d) = %f \n",i+1,getinput4_.wavobs[i],i,fwhm[i]);
795  }
796
797  P->idx450 = bindex_get(450); // get the index nearest the 450nm wavelength
798  if (P->idx450 < 0) {
799  printf("-E- %s line %d : Unable to determine 450 nm index from spectrum.\n",
800  __FILE__, __LINE__);
801  exit(FATAL_ERROR);
802  }
803  // for (i=0;i<nbands;i++) {
804  // getinput4_.fwhm[i] = fwhm[i];
805  // // printf("->RJH:fwave(%d)=%f fwhm(%d) = %f\n",i+1,getinput4_.wavobs[i],i,fwhm[i]);
806  // }
807  /*
808  // Determine which atmospheric model to use
809  // call to get_input sets h,p,t, and vmr arrays to default model 1
810  // 1 = tropical
811  // 2 = mid latitude summer
812  // 3 = mid latitude winter
813  // 4 = subarctic summer
814  // 5 = subarctic winter
815  // 6 = US standard 1962
816  // 7 = user defined model
817
818  //INIT
819  //Get Atmospheric model to use
820  */
821  *model = 6;
822
823  if (atrem_model == 0) {
824  int16_t year, day;
825  double secs;
826  unix2yds(l1rec->scantime, &year, &day, &secs);
827  *model = getModelNum(*l1rec->lat, day);
828  } else {
829  if (atrem_model <= 6) {
830  *model = atrem_model;
831  // nwave = rdatreminfo(sensorID, 0, "model", (void **)&model);
832  } else {
833  printf("-E- %s line %d : Invalid atmospheric model, atrem_model = %d. Valid range is 0-6.\n",
834  __FILE__, __LINE__, atrem_model);
835  exit(FATAL_ERROR);
836  }
837  }
838  P->model = *model;
839  nb = init_tpvmr(P->model);
840
841  if (nb <= 0) {
842  printf("-E- %s line %d : Atmospheric data could not be initialized.\n",
843  __FILE__, __LINE__);
844  exit(FATAL_ERROR);
845  }
846
847  getinput5_.nobs = nbands; //l1rec->nbands;
848  getinput5_.dlt = 0; //used internally to the fortran routines
849  getinput5_.dlt2 = 0; //to determine what to use for the width of the bands
850  getinput5_.hsurf = 0; //mean surface elevation
851  getinput14_.xppp = 700; //Plane/Satellite altitude (km)
852
853  if (getinput5_.hsurf < 0 || getinput5_.hsurf > getinput3_.h[nb - 2]) {
854  printf("-E- %s line %d : Elevation=%f must be less than the maximum elevation in the atmospheric model.\n",
855  __FILE__, __LINE__, getinput5_.hsurf);
856  exit(FATAL_ERROR);
857  }
858
859  nwave = rdatreminfo(sensorID, 0, "xppp", (void **) &xppp);
860
861  if (getinput14_.xppp < getinput5_.hsurf) {
862  printf("-E- %s line %d : Sensor altitude=%f must be greater than the bottom surface elevation.\n",
863  __FILE__, __LINE__, getinput14_.xppp);
864  exit(FATAL_ERROR);
865
866  }
867  /*
868  // Get ranges on curve for the two atmospheric windows surrounding the 1.14-um
869  // water vapor absorption feature and the center point of the 1.14-um water
870  // vapor absorption feature. Enter:
871  // 1. the midpoint of third window (0.6-2.5)
872  // 2. number of points to average for third window (1-10)
873  // 3. the midpoint of fourth window (0.6-2.5)
874  // 4. number of points to average for fourth window (1-10)
875  // 5. the midpoint of 1.14-um absorption feature (0.6-2.5)
876  // 6. the number of points to average for the absorption feature (1-30)
877  */
878  // This is the default for HICO HS
879  float *wndow1 = 0, *wndow2 = 0, *wp94c = 0, *wndow3 = 0, *wndow4 = 0, *w1p14c = 0;
880
881  wndow1 = (float *) calloc(1, sizeof (float));
882  wndow2 = (float *) calloc(1, sizeof (float));
883  wndow3 = (float *) calloc(1, sizeof (float));
884  wndow4 = (float *) calloc(1, sizeof (float));
885  wp94c = (float *) calloc(1, sizeof (float));
886  w1p14c = (float *) calloc(1, sizeof (float));
887
888  getinput6_.wndow1 = 0.705;
889  getinput6_.wndow2 = 0.745;
890  getinput6_.wp94c = 0.725;
891  getinput6_.wndow3 = 0.805;
892  getinput6_.wndow4 = 0.845;
893  getinput6_.w1p14c = 0.825;
894
895  nwave = rdatreminfo(sensorID, 0, "window1", (void **) &wndow1);
896  nwave = rdatreminfo(sensorID, 0, "window2", (void **) &wndow2);
897  nwave = rdatreminfo(sensorID, 0, "window3", (void **) &wndow3);
898  nwave = rdatreminfo(sensorID, 0, "window4", (void **) &wndow4);
899  nwave = rdatreminfo(sensorID, 0, "wp94c", (void **) &wp94c);
900  nwave = rdatreminfo(sensorID, 0, "w1p14c", (void **) &w1p14c);
901
902  getinput6_.wndow1 = *wndow1;
903  getinput6_.wndow2 = *wndow2;
904  getinput6_.wp94c = *wp94c;
905  getinput6_.wndow3 = *wndow3;
906  getinput6_.wndow4 = *wndow4;
907  getinput6_.w1p14c = *w1p14c;
908
909  if (getinput6_.wndow1 < 0.6 || getinput6_.wndow1 > 2.5) {
910  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
911  "in the .94-um region1:%f\nValid values are 0.6-2.5\n", getinput6_.wndow1);
912  printf("-E- %s line %d : See stderr output for details.\n",
913  __FILE__, __LINE__);
914  exit(FATAL_ERROR);
915  }
916  if (getinput6_.wndow2 < 0.6 || getinput6_.wndow2 > 2.5) {
917  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
918  "in the .94-um region2:%f\nValid values are 0.6-2.5\n", getinput6_.wndow2);
919  printf("-E- %s line %d : See stderr output for details.\n",
920  __FILE__, __LINE__);
921  exit(FATAL_ERROR);
922  }
923  if (getinput6_.wp94c <= getinput6_.wndow1 || getinput6_.wp94c >= getinput6_.wndow2) {
924  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
925  "in the .94-um region:%f\nValid range is: %f < value < %f \n", getinput6_.wp94c, getinput6_.wndow1, getinput6_.wndow2);
926  printf("-E- %s line %d : See stderr output for details.\n",
927  __FILE__, __LINE__);
928  exit(FATAL_ERROR);
929  }
930  if (getinput6_.wndow3 < 0.6 || getinput6_.wndow3 > 2.5) {
931  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
932  "in the .94-um region3:%f\nValid values are 0.6-2.5\n", getinput6_.wndow3);
933  printf("-E- %s line %d : See stderr output for details.\n",
934  __FILE__, __LINE__);
935  exit(FATAL_ERROR);
936  }
937  if (getinput6_.wndow4 < 0.6 || getinput6_.wndow4 > 2.5) {
938  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
939  "in the .94-um region4:%f\nValid values are 0.6-2.5\n", getinput6_.wndow4);
940  printf("-E- %s line %d : See stderr output for details.\n",
941  __FILE__, __LINE__);
942  exit(FATAL_ERROR);
943  }
944  if (getinput6_.w1p14c <= getinput6_.wndow3 || getinput6_.w1p14c >= getinput6_.wndow4) {
945  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
946  "in the .94-um region:%f\nValid range is: %f < value < %f \n", getinput6_.w1p14c, getinput6_.wndow3, getinput6_.wndow4);
947  printf("-E- %s line %d : See stderr output for details.\n",
948  __FILE__, __LINE__);
949  exit(FATAL_ERROR);
950  }
951
952  // Set full_calc to 1 to turn on explicit, full calculation of
953  // water vapor correction. Setting value to 0 turns on speedy
954  // h2o absorption calculation and sets all other gas options to 0
955
956  // int32_t *full_calc;
957  // full_calc = (int32_t *) calloc(1,sizeof(int32_t));
958  // if (atrem_full != 0) {
959  // *full_calc = (atrem_full < 0? 0:1);
960  // } else {
961  // nwave = rdatreminfo(sensorID, 0, "full_calc", (void **) &full_calc);
962  // }
963
964  if (atrem_splitpaths > 0 && atrem_full == 0) {
965  printf("ATREM: Turning on atrem_full because you selected atrem_splitpaths\n");
966  atrem_full = 1;
967  }
968  if (atrem_full > 0)
969  printf("ATREM: Warning : full_calc !=0. Atrem will calculate transmittance table for every pixel\n");
970
971  // getinput5_.full_calc = *full_calc;
972
973  getinput5_.full_calc = atrem_full;
974
975  // int32_t *dogeom;
976  // dogeom = (int32_t *) calloc(1,sizeof(int32_t));
977  // if (atrem_geom != 0) {
978  // *dogeom = (atrem_geom < 0? 0:1);
979  // } else {
980  // nwave = rdatreminfo(sensorID, 0, "dogeom", (void **) &dogeom);
981  // }
982  //
983  // P->dogeom = *dogeom;
984
985  if (atrem_geom > 0)
986  printf("ATREM: Warning : dogeom !=0. Geometry will be calculated every pixel\n");
987
988  P->dogeom = atrem_geom;
989  //Default for HICO HS
990  int32_t *nb1, *nb2, *nb3, *nb4, *nbp94, *nb1p14;
991  nb1 = (int32_t *) calloc(1, sizeof (int32_t));
992  nb2 = (int32_t *) calloc(1, sizeof (int32_t));
993  nb3 = (int32_t *) calloc(1, sizeof (int32_t));
994  nb4 = (int32_t *) calloc(1, sizeof (int32_t));
995  nbp94 = (int32_t *) calloc(1, sizeof (int32_t));
996  nb1p14 = (int32_t *) calloc(1, sizeof (int32_t));
997
998  getinput7_.nb1 = 3;
999  getinput7_.nb2 = 3;
1000  getinput7_.nb3 = 3;
1001  getinput7_.nb4 = 3;
1002  getinput7_.nbp94 = 5;
1003  getinput7_.nb1p14 = 5;
1004
1005  nwave = rdatreminfo(sensorID, 0, "nb1", (void **) &nb1);
1006  nwave = rdatreminfo(sensorID, 0, "nb2", (void **) &nb2);
1007  nwave = rdatreminfo(sensorID, 0, "nb3", (void **) &nb3);
1008  nwave = rdatreminfo(sensorID, 0, "nb4", (void **) &nb4);
1009  nwave = rdatreminfo(sensorID, 0, "nbp94", (void **) &nbp94);
1010  nwave = rdatreminfo(sensorID, 0, "nb1p14", (void **) &nb1p14);
1011
1012  getinput7_.nb1 = *nb1;
1013  getinput7_.nb2 = *nb2;
1014  getinput7_.nb3 = *nb3;
1015  getinput7_.nb4 = *nb4;
1016  getinput7_.nbp94 = *nbp94;
1017  getinput7_.nb1p14 = *nb1p14;
1018
1019  if (getinput7_.nb1 < 1 || getinput7_.nb1 > 50) {
1020  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1021  "the .94-um region:%d\nValid values are 1-50\n", getinput7_.nb1);
1022  printf("-E- %s line %d : See stderr output for details.\n",
1023  __FILE__, __LINE__);
1024  exit(FATAL_ERROR);
1025  }
1026
1027  if (getinput7_.nb2 < 1 || getinput7_.nb2 > 50) {
1028  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1029  "the .94-um region:%d\nValid values are 1-50\n", getinput7_.nb2);
1030  printf("-E- %s line %d : See stderr output for details.\n",
1031  __FILE__, __LINE__);
1032  exit(FATAL_ERROR);
1033  }
1034  if (getinput7_.nbp94 < 1 || getinput7_.nbp94 > 90) {
1035  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1036  "the .94-um region:%d\nValid values are 1-90\n", getinput7_.nbp94);
1037  printf("-E- %s line %d : See stderr output for details.\n",
1038  __FILE__, __LINE__);
1039  exit(FATAL_ERROR);
1040  }
1041  if (getinput7_.nb3 < 1 || getinput7_.nb3 > 50) {
1042  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1043  "the .94-um region:%d\nValid values are 1-50\n", getinput7_.nb3);
1044  printf("-E- %s line %d : See stderr output for details.\n",
1045  __FILE__, __LINE__);
1046  exit(FATAL_ERROR);
1047  }
1048  if (getinput7_.nb4 < 1 || getinput7_.nb4 > 50) {
1049  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1050  "the .94-um region:%d\nValid values are 1-50\n", getinput7_.nb4);
1051  printf("-E- %s line %d : See stderr output for details.\n",
1052  __FILE__, __LINE__);
1053  exit(FATAL_ERROR);
1054  }
1055  if (getinput7_.nb1p14 < 1 || getinput7_.nb1p14 > 110) {
1056  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1057  "the .94-um region:%d\nValid values are 1-90\n", getinput7_.nbp94);
1058  printf("-E- %s line %d : See stderr output for details.\n",
1059  __FILE__, __LINE__);
1060  exit(FATAL_ERROR);
1061  }
1062
1063  //INIT
1064  int32_t *h2o, *co2, *o3, *n2o, *co, *ch4, *o2, *no2;
1065  h2o = (int32_t *) calloc(1, sizeof (int32_t));
1066  co2 = (int32_t *) calloc(1, sizeof (int32_t));
1067  o3 = (int32_t *) calloc(1, sizeof (int32_t));
1068  n2o = (int32_t *) calloc(1, sizeof (int32_t));
1069  co = (int32_t *) calloc(1, sizeof (int32_t));
1070  ch4 = (int32_t *) calloc(1, sizeof (int32_t));
1071  o2 = (int32_t *) calloc(1, sizeof (int32_t));
1072  no2 = (int32_t *) calloc(1, sizeof (int32_t));
1073
1074  *h2o = 1;
1075
1076  getinput1_.h2o = 1;
1077  getinput1_.co2 = 0;
1078  getinput1_.o3 = 0;
1079  getinput1_.n2o = 0;
1080  getinput1_.co = 0;
1081  getinput1_.ch4 = 0;
1082  getinput1_.o2 = 0;
1083  getinput1_.no2 = 0;
1084
1085  // nwave = rdatreminfo(sensorID, 0, "h2o", (void **) &h2o);
1086  nwave = rdatreminfo(sensorID, 0, "co2", (void **) &co2);
1087  nwave = rdatreminfo(sensorID, 0, "o3", (void **) &o3);
1088  nwave = rdatreminfo(sensorID, 0, "n2o", (void **) &n2o);
1089  nwave = rdatreminfo(sensorID, 0, "co", (void **) &co);
1090  nwave = rdatreminfo(sensorID, 0, "ch4", (void **) &ch4);
1091  nwave = rdatreminfo(sensorID, 0, "o2", (void **) &o2);
1092  nwave = rdatreminfo(sensorID, 0, "no2", (void **) &no2);
1093
1094
1095  getinput1_.co2 = *co2;
1096  getinput1_.o3 = *o3;
1097  getinput1_.n2o = *n2o;
1098  getinput1_.co = *co;
1099  getinput1_.ch4 = *ch4;
1100  getinput1_.o2 = *o2;
1101  getinput1_.no2 = *no2;
1102
1103  // The typical ozone amount is: 0.28-0.55 (atm-cm). The built-in NO2
1104  // column amount is 5.0E+15 molecules/cm^2.
1105  if (atrem_opt > 0) {
1106  getinput1_.co2 = (atrem_opt & ATREM_CO2) > 0;
1107  getinput1_.o3 = (atrem_opt & ATREM_O3) > 0;
1108  getinput1_.n2o = (atrem_opt & ATREM_N2O) > 0;
1109  getinput1_.co = (atrem_opt & ATREM_CO) > 0;
1110  getinput1_.ch4 = (atrem_opt & ATREM_CH4) > 0;
1111  getinput1_.o2 = (atrem_opt & ATREM_O2) > 0;
1112  getinput1_.no2 = (atrem_opt & ATREM_NO2) > 0;
1113
1114  }
1115
1116  printf("ATREM Gas Options:H2O:1 CO2:%d O3:%d N2O:%d CO:%d CH4:%d O2:%d NO2:%d\n",
1117  getinput1_.co2, getinput1_.o3, getinput1_.n2o, getinput1_.co, getinput1_.ch4, getinput1_.o2, getinput1_.no2);
1118
1119  // Now check to make sure the same gas option isn't selected from the l2gen gas option.
1120
1121  flag_gas = 0;
1122  if (gas_opt & H2O_BIT) {
1123  printf("ATREM: cannot be used with gas_opt bit mask=%d (H2O)\n", H2O_BIT);
1124  flag_gas = 1;
1125  }
1126  if (getinput1_.no2 == 1 && (gas_opt & NO2_BIT)) {
1127  printf("ATREM: cannot be used with gas_opt bit mask=%d (NO2)\n", NO2_BIT);
1128  flag_gas = 1;
1129
1130  }
1131  if (getinput1_.co2 == 1 && (gas_opt & CO2_BIT)) {
1132  printf("ATREM: cannot be used with gas_opt bit mask=%d (CO2)\n", CO2_BIT);
1133  flag_gas = 1;
1134  }
1135  if (getinput1_.o3 == 1 && (gas_opt & O3_BIT)) {
1136  printf("ATREM: cannot be used with gas_opt bit mask=%d (O3)\n", O3_BIT);
1137  flag_gas = 1;
1138  }
1139
1140  if (flag_gas) {
1141  printf("Error: Conflict using ATREM (gas_opt=16 bitmask) with atrem_opt=%d and gas_opt=%d. " \
1142  "\nPlease resolve. Refer to command line options for atrem_opt and gas_opt.\n", atrem_opt, gas_opt);
1143  exit(1);
1144  }
1145
1146  float *vrto3, *sno2;
1147  vrto3 = (float *) calloc(1, sizeof (float));
1148  sno2 = (float *) calloc(1, sizeof (float));
1149
1150  getinput3_.vrto3 = 0.34; //total column ozone amount (atm-cm)
1151  getinput3_.sno2 = 1.0; //NO2 scaling factor (to 5.E+15 molecules/cm^2)
1152
1153  nwave = rdatreminfo(sensorID, 0, "vrto3", (void **) &vrto3);
1154  nwave = rdatreminfo(sensorID, 0, "sno2", (void **) &sno2);
1155
1156  getinput3_.vrto3 = *vrto3; //total column ozone amount (atm-cm)
1157  getinput3_.sno2 = *sno2; //NO2 scaling factor (to 5.E+15 molecules/cm^2)
1158
1160
1161  P->nb1 = getinput7_.nb1;
1162  P->nb2 = getinput7_.nb2;
1163  P->nb3 = getinput7_.nb3;
1164  P->nb4 = getinput7_.nb4;
1165  P->nbp94 = getinput7_.nbp94;
1166  P->nb1p14 = getinput7_.nb1p14;
1167  P->nobs = getinput5_.nobs;
1168  P->delta = getinput5_.dlt;
1169  P->delta2 = getinput5_.dlt2;
1170
1171  return nwave;
1172 }
1173
1174 int get_angle_limits(float **anglelimit, float **insenz, float **insolz, int *n_senz, int *n_solz) {
1175
1176  char filename[FILENAME_MAX];
1177  char *infile = "atrem_angle_limit.nc";
1178  char *filedir;
1179
1180  char name[H4_MAX_NC_NAME];
1181  char sdsname[H4_MAX_NC_NAME];
1182  int ncid, ndims;
1183  int32 sds_id;
1184  int status;
1185  nc_type rh_type; /* variable type */
1186  int dimids[H4_MAX_VAR_DIMS]; /* dimension IDs */
1187  int natts; /* number of attributes */
1188  size_t length;
1189
1190  static float *senz, *solz;
1191  static float *angle_limit;
1192
1193  if ((filedir = getenv("OCDATAROOT")) == NULL) {
1194  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
1195  return (-1);
1196  }
1197  strcpy(filename, filedir);
1198  strcat(filename, "/common/");
1199  strcat(filename, infile);
1200  strcat(filename, "\0");
1201
1202  printf("ATREM Angle_Limit FILE=%s\n", filename);
1203
1204  if (nc_open(filename, NC_NOWRITE, &ncid) == NC_NOERR) {
1205
1206  strcpy(sdsname, "senz");
1207
1208  status = nc_inq_varid(ncid, sdsname, &sds_id);
1209  if (status != NC_NOERR) {
1210  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1211  __FILE__, __LINE__, sdsname, infile);
1212  exit(1);
1213  }
1214
1215  status = nc_inq_var(ncid, sds_id, 0, &rh_type, &ndims, dimids,
1216  &natts);
1217
1218  if (nc_inq_dimlen(ncid, dimids[0], &length) != NC_NOERR) {
1219  nc_inq_dim(ncid, dimids[0], name, &length);
1220  fprintf(stderr,
1221  "-E- %s line %d: could not get size of dimension \"%s\" in netCDF File.\n",
1222  __FILE__, __LINE__, name);
1223  exit(1);
1224  }
1225
1226  *n_senz = length;
1227
1228  if ((senz = (float *) calloc(*n_senz, sizeof (float))) == NULL) {
1229  printf("-E- : Error allocating memory to senz\n");
1230  exit(FATAL_ERROR);
1231  }
1232
1233  if (nc_get_var(ncid, sds_id, senz) != NC_NOERR) {
1234  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1235  __FILE__, __LINE__, sdsname, infile);
1236  exit(1);
1237  }
1238
1239  strcpy(sdsname, "solz");
1240
1241  status = nc_inq_varid(ncid, sdsname, &sds_id);
1242  if (status != NC_NOERR) {
1243  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1244  __FILE__, __LINE__, sdsname, infile);
1245  exit(1);
1246  }
1247
1248  status = nc_inq_var(ncid, sds_id, 0, &rh_type, &ndims, dimids,
1249  &natts);
1250  if (status != NC_NOERR) {
1251  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1252  __FILE__, __LINE__, sdsname, infile);
1253  exit(1);
1254  }
1255
1256  if (nc_inq_dimlen(ncid, dimids[0], &length) != NC_NOERR) {
1257  nc_inq_dim(ncid, dimids[0], name, &length);
1258  fprintf(stderr,
1259  "-E- %s line %d: could not get size of dimension \"%s\" in netCDF File.\n",
1260  __FILE__, __LINE__, name);
1261  exit(1);
1262  }
1263
1264  *n_solz = length;
1265
1266  if ((solz = (float *) calloc(*n_solz, sizeof (float))) == NULL) {
1267  printf("-E- : Error allocating memory to solz\n");
1268  exit(FATAL_ERROR);
1269  }
1270
1271  if (nc_get_var(ncid, sds_id, solz) != NC_NOERR) {
1272  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1273  __FILE__, __LINE__, sdsname, infile);
1274  exit(1);
1275  }
1276
1277  strcpy(sdsname, "angle_limit");
1278
1279  status = nc_inq_varid(ncid, sdsname, &sds_id);
1280  if (status != NC_NOERR) {
1281  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1282  __FILE__, __LINE__, sdsname, infile);
1283  exit(1);
1284  }
1285
1286  status = nc_inq_var(ncid, sds_id, 0, &rh_type, &ndims, dimids,
1287  &natts);
1288  if (status != NC_NOERR) {
1289  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1290  __FILE__, __LINE__, sdsname, infile);
1291  exit(1);
1292  }
1293  if (ndims != 2) {
1294  fprintf(stderr, "-E- %s line %d: Wrong number of dimensions for %s. Need 2 got %d.\n",
1295  __FILE__, __LINE__, sdsname, ndims);
1296  exit(1);
1297
1298  }
1299  // printf("ANGLE_LIMIT: %s n_senz=%d n_solz=%d\n",infile,*n_senz,*n_solz);
1300
1301  if ((angle_limit = (float *) calloc((*n_senz) * (*n_solz), sizeof (float))) == NULL) {
1302  printf("-E- : Error allocating memory to angle_limit\n");
1303  exit(FATAL_ERROR);
1304  }
1305
1306  // for (i=0; i< *n_senz;i++) {
1307  // if ( (anglelimit[i] = (float *)calloc((*n_solz) ,sizeof(float))) == NULL) {
1308  // printf("-E- : Error allocating memory to tindx\n");
1309  // exit(FATAL_ERROR);
1310  // }
1311  // }
1312
1313  if (nc_get_var(ncid, sds_id, angle_limit) != NC_NOERR) {
1314  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1315  __FILE__, __LINE__, sdsname, infile);
1316  exit(1);
1317  }
1318
1319  // for (j=0;j<*n_senz;j++)
1320  // for (i=0;i<*n_solz;i++) {
1321  // //anglelimit[j][i] = *(angle_limit+j*(*n_solz)+i);
1322  // printf("RJH: angle_limt: %f %f %f\n",solz[i],senz[j],*(angle_limit+j*(*n_solz)+i));
1323  //
1324  // }
1325
1326  *insenz = (float *) senz;
1327  *insolz = (float *) solz;
1328  *insolz = (float *) solz;
1329  *anglelimit = (float *) angle_limit;
1330
1331  } else {
1332  fprintf(stderr, "-E- %s line %d: Error opening infile = %s.\n",
1333  __FILE__, __LINE__, infile);
1334  return (1);
1335  }
1336
1337
1338  return (0);
1339
1340 }
1341
1342 float get_current_angle_limit(float insenz, float insolz, int *ii, int *jj, float **anglelimit, float *senz, float *solz, int n_senz, int n_solz) {
1343  int i, j;
1344
1345  i = *ii;
1346  j = *jj;
1347
1348  if (i < 0 || i >= n_senz) i = 0;
1349  if (j < 0 || j >= n_solz) j = 0;
1350
1351  if (insenz > senz[i]) {
1352  while (i < n_senz && insenz > senz[i])
1353  i++;
1354  if (i >= n_senz) i = n_senz - 1;
1355  } else {
1356  while (i >= 0 && insenz < senz[i])
1357  i--;
1358  if (i < 0) i = 0;
1359  }
1360
1361  if (insolz > solz[j]) {
1362  while (j < n_solz && insolz > solz[j])
1363  j++;
1364  if (j >= n_solz) j = n_solz - 1;
1365  } else {
1366  while (j >= 0 && insolz < solz[j])
1367  j--;
1368  if (j < 0) j = 0;
1369  }
1370
1371  *ii = i;
1372  *jj = j;
1373
1374  return (anglelimit[i][j]);
1375 }
1376
1377 /********************************************************************************
1378 !* *
1379 !* Name: FINDMATCH *
1380 !* Purpose: finds the closest match for ELEM in LIST *
1381 !* Parameters: LIST - array of values to match ELEM to. Elements is array *
1382 !* should increase in value. *
1383 !* Algorithm: linearly compare ELEM to each element. When ELEM is smaller *
1384 !* than the LIST(I), then it is assumed to be closest to LIST(I-1) *
1385 !* or LIST(I). The one that has the smallest absolute difference *
1386 !* to ELEM is returned as the closest match. *
1387 !* Globals used: none *
1388 !* Global output: none *
1389 !* Return Codes: the closest matching element index *
1390 !* Special Considerations: none *
1391 !* *
1392 !********************************************************************************/
1393
1394 int findMatch_(float *list, int *nobs, float *elem) {
1395  int i;
1396  float diff1, diff2;
1397
1398  for (i = 0; i<*nobs && list[i] <= *elem; i++)
1399
1400  diff1 = fabs(list[i - 1] - *elem);
1401  diff2 = fabs(list[i] - *elem);
1402
1403  if (diff1 < diff2)
1404  return (i - 1);
1405  else
1406  return (i);
1407
1408 }
#define ATREM_N2O
Definition: atrem_corl1.h:28
float tg_sen[NBANDS]
Definition: atrem_corl1.h:127
integer, parameter int16
Definition: cubeio.f90:3
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:217
int32_t co
Definition: atrem_corl1.h:118
struct @25 getinput4_
#define H2O_BIT
Definition: l12_parms.h:96
int j
Definition: decode_rs.h:73
void init_speccal_()
int32_t day
int status
Definition: l1_czcs_hdf.c:31
struct @24 getinput3_
struct @22 input_l2gen_
int32_t o3
Definition: atrem_corl1.h:118
float wndow4
Definition: atrem_corl1.h:153
int32_t co2
Definition: atrem_corl1.h:118
#define NULL
Definition: decode_rs.h:63
float sno2
Definition: atrem_corl1.h:133
struct @9 tran_table1_
float wndow1
Definition: atrem_corl1.h:153
float wndow2
Definition: atrem_corl1.h:153
#define CO2_BIT
Definition: l12_parms.h:94
void geometry_()
int32_t nb4
Definition: atrem_cor.h:98
list(APPEND LIBS \${NETCDF_LIBRARIES}) list(APPEND LIBS \$
Definition: CMakeLists.txt:9
float * lat
#define O3_BIT
Definition: l12_parms.h:93
int32_t nbp94
Definition: atrem_cor.h:98
int init_atrem(int32_t sensorID, paramstr *P, l1str *l1rec, int32_t nbands)
struct @46 debug_atrem
int bindex_get(int32_t wave)
Definition: windex.c:43
float fwhm[NBANDS]
Definition: atrem_corl1.h:144
struct @30 getinput14_
struct @4 init_speccal6_
float vrto3
Definition: atrem_corl1.h:133
int32_t nobs
Definition: atrem_cor.h:93
instr * input
#define MODELMAX
Definition: atrem_corl1.h:20
int32_t h2o
Definition: atrem_corl1.h:118
int jb
Definition: atrem_corl1.h:227
struct @6 init_speccal8_
int findMatch_(float *list, int *nobs, float *elem)
int32_t no2
Definition: atrem_corl1.h:118
#define ATREM_NO2
Definition: atrem_corl1.h:24
#define NH2OMAX
Definition: atrem_corl1.h:17
struct @21 getinput1_
int32_t nb1
Definition: atrem_cor.h:98
struct @7 init_speccal10_
#define ATREM_O3
Definition: atrem_corl1.h:22
int32_t nb2
Definition: atrem_cor.h:98
float get_current_angle_limit(float insenz, float insolz, int *ii, int *jj, float **anglelimit, float *senz, float *solz, int n_senz, int n_solz)
#define ATREM_CO
Definition: atrem_corl1.h:25
struct @1 getinput5_
int getModelNum(float lat, int32_t day)
float tg_sol[NBANDS]
Definition: atrem_corl1.h:127
#define NO2_BIT
Definition: l12_parms.h:95
#define FATAL_ERROR
Definition: swl0_parms.h:5
float w1p14c
Definition: atrem_corl1.h:153
const double delta
void unix2yds(double usec, short *year, short *day, double *secs)
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
int32_t o2
Definition: atrem_corl1.h:118
float ja
Definition: atrem_cor.h:114
#define ATREM_CH4
Definition: atrem_corl1.h:26
int init_tpvmr(int model)
float wp94c
Definition: atrem_corl1.h:153
void get_input_()
int32_t hunt(float *xx, int32_t n, double x, int32_t jlo)
#define ATREM_O2
Definition: atrem_corl1.h:27
float get_atrem(float *tg_tot, float *rhot, paramstr *P)
int get_atrem_cor(int32_t sensorID, l1str *l1rec, int32_t ip, float *rhot, float *tg_tot, float *tg_sol, float *tg_sen)
struct @44 geometry_l2gen_
int32_t nbands
@ NBANDS
Definition: make_L3_v1.1.c:53
struct @27 getinput6_
#define fabs(a)
Definition: misc.h:93
char * name
Definition: Granule.c:1234
float xppp
Definition: atrem_corl1.h:165
struct @10 geometry3_
int32_t ch4
Definition: atrem_corl1.h:118
struct @45 tpvmr_init1_
struct @23 tran_table_l2gen_
float wndow3
Definition: atrem_corl1.h:153
int32_t rdatreminfo(int32_t sensorID, int32_t evalmask, const char *pname, void **pval)
Definition: rdatreminfo.c:38
int32_t model
Definition: atrem_corl1.h:132
struct @3 init_speccal3_
#define ATREM_CO2
Definition: atrem_corl1.h:23
int get_angle_limits(float **anglelimit, float **insenz, float **insolz, int *n_senz, int *n_solz)
int i
Definition: decode_rs.h:71
struct @8 init_speccal11_
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:95
int32_t nb1p14
Definition: atrem_cor.h:98
struct @5 init_speccal7_
#define NH2OMAXM1
Definition: atrem_corl1.h:18
int32_t nb3
Definition: atrem_cor.h:98
struct @2 getinput7_
int32_t nb
Definition: atrem_corl1.h:132
int32_t n2o
Definition: atrem_corl1.h:118
void tran_table_()