OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
get_atrem_corl1v3.c
Go to the documentation of this file.
1 /*
2  * get_atrem_corl1v3.c
3  *
4  * ATmospheric REMoval (ATREM)
5  *
6  * Created on: Feb 19, 2015
7  * Author: Rick Healy (richard.healy@nasa.gov)
8  * from ATREM Fortran Code by Bo Cai Gao
9  * Comments from Fortran Code unless noted otherwise
10  * Variable names may differ slightly from comments
11  * (mostly changed to lower case)
12  * nbands = NOBS in Fortran code
13  *
14  * 05/05/2016 - C stand alone (Version 3) completed for the most part - clean up continues...
15  * 03/01/2016 - Began converting rest of code to C
16  * Netcdf files created
17  * 01/15/2016 - Split the solar and sensor paths directly inside atrem.
18  * specify with atrem_splitpaths. Only works with atrem_full on.
19  *
20  * 10/20/2015 - r.healy - added ATREM options to command line
21  * atrem_opt - select gases for transmittance calculation
22  * atrem_geom - turn on/off geometry recalculation
23  * 0 = geometry calculated on error angle limit
24  * 1 = recalculate every pixel
25  * atrem_full - turn on/off full calculation (k-dist otherwise)
26  * atrem_model - select atmospheric model (1-6)
27  * - 0= determine from latitude and day of year
28  *
29  * * Notes about water vapor VMRS and related quantities: *
30  * *
31  * VAPVRT(60) - a table containing 60 column vapor values (in unit of cm) *
32  * *
33  * VAP_SLANT(I) = VAPVRT(I) * 2.0, VAP_SLANT is a new table for containing *
34  * two-way total vapor amounts. Here the number "2" can be *
35  * changed to other numbers, e.g., 2.5, without major *
36  * effects on retrieved water vapor values. *
37  * *
38  * G_VAP(I = 1,..., NL) = true vapor geometric factor for each layer in *
39  * the model atmosphere (after adjusting for the elevated *
40  * surface. *
41  * *
42  * VMRM(I) = VMRM(I)*G_VAP(I). The VMRS are multiplied by the geometrical *
43  * factor. We can calculate the vapor transmittance on the *
44  * Sun-surface-sensor path by assuming a vertical path in *
45  * the model atmosphere with geometric-factor-adjusted VMRS. *
46  * *
47  * CLMVAP = vertical column amount from ground to space in model atmosphere*
48  * CLMVAPP = vertical column amount from ground to aircraft or satellite *
49  * sensor in model atmosphere *
50  * Q = 2.152E25 = # of molecules above the surface at one atmosphere *
51  * (in unit of molecules/cm**2) *
52  * *
53  * VAP_SLANT_MDL= CLMVAP/COS(SOLZNI) + CLMVAPP/COS(OBSZNI) = total amount *
54  * of water vapor in the model atmosphere in the L-shaped *
55  * Sun-surface-plane ray path. *
56  * *
57  * G_VAP_EQUIV = VAP_SLANT_MDL / CLMVAP = the "equivalent" geometrical *
58  * factor corresponding to the total slant vapor amount *
59  * VAP_SLANT_MDL and the column vapor amount CLMVAP. *
60  * *
61  * SSH2O(I) (I = 1, ..., 60) - a pure scaling factor relative to the total *
62  * slant vapor amount of VAP_SLANT_MDL, and *
63  * SSH2O(I) = VAP_SLANT(I) / VAP_SLANT_MDL *
64  * *
65  * SH2O = one value of SSH2O(I). SH2O is used during generation of the *
66  * look-up table. *
67  * *
68  * VAPTT = VAP_SLANT_MDL*SH2O, is the absolute total vapor amount on the *
69  * L-shaped path corresponding to a spectrum stored in the *
70  * look-up table. *
71  * *
72  * CLMWVP = 0.5*(VAPTTA+VAPTTB)/G_VAP_EQUIV, is the retrieved column water *
73  * vapor amount from imaging spectrometer data. *
74  ********************************************************************************
75  ********************************************************************************
76  ********************************************************************************!*
77 References: ********************************************************************!* *
78 !* Gao, B.-C., K. Heidebrecht, and A. F. H. Goetz, Derivation of scaled *
79 !* surface reflectances from AVIRIS data, Remote Sens. Env., 44, *
80 !* 165-178, 1993. *
81 !* Gao, B.-C., and C. O. Davis, Development of an operational algorithm *
82 !* for removing atmospheric effects from HYDICE and HSI data, *
83 !* in SPIE'96 Conference Proceedings, Vol. 2819, 45-55, 1996. *
84 !* Gao, B.-C., and A. F. H. Goetz, Column atmospheric water vapor and *
85 !* vegetation liquid water retrievals from airborne imaging *
86 !* spectrometer data, J. Geophys. Res., 95, 3549-3564, 1990. *
87 !* Goetz, A. F. H., and M. Herring, The high resolution imaging *
88 !* spectrometer (HIRIS) for Eos, IEEE Trans. Geosci. Remote Sens., 27,*
89 !* 136-144, 1989. *
90 !* Goetz, A. F. H., G. Vane, J. Solomon, and B. N. Rock, Imaging *
91 !* spectrometry for Earth remote sensing, Science, 228, 1147-1153,1985*
92 !* Green, R. O., and B.-C. Gao, A Proposed Update to the Solar Irradiance*
93 !* Spectrum used in LOWTRAN and MODTRAN, in Summaries of the Fourth *
94 !* Annual JPL Airborne Geoscience Workshop, October 25-29, (Editor, *
95 !* R. O. Green), JPL Publ. 93-26, Vol. 1, pp. 81-84, Jet Propul. Lab, *
96 !* Pasadena, Calif., 1993. *
97 !* Kneizys, F. X., E. P. Shettle, L. W. Abreu, J. H. Chetwynd, G. P. *
98 !* Anderson, W. O. Gallery, J. E. A. Selby, and S. A. Clough, Users *
99 !* guide to LOWTRAN7, AFGL-TR-8-0177, Air Force Geophys. Lab., *
100 !* Bedford, Mass., 1988. *
101 !* Iqbal, M., An Introduction To Solar Radiation, Academic, San Diego, *
102 !* Calif., 1983. *
103 !* Malkmus, W., Random Lorentz band model with exponential-tailed S line *
104 !* intensity distribution function, J. Opt. Soc. Am., 57, 323-329,1967*
105 !* Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, *
106 !* Numerical Recipes-The ART of Scientific Computing, Cambridge *
107 !* University Press, 1986. *
108 !* Rothman, L. S., et al., The HITRAN 2008 molecular spectroscopic *
109 !* database, JQSRT, 110, 533-572, 2009. *
110 !* Solomon, S., R. W. Portmann, R. W. Sanders, J. S. Daniel, W. Madsen, *
111 !* B. Bartram, and E. G. Dutton, On the role of nitrogen dioxide in *
112 !* the absorption of solar radiation, J. Geophys. Res., 104, *
113 !* 12047-12058, 1999. *
114 !* Tanre, D., C. Deroo, P. Duhaut, M. Herman, J. J. Morcrette, J. Perbos,*
115 !* and P. Y. Deschamps, Description of a computer code to simulate *
116 !* the satellite signal in the solar spectrum: the 5S code, Int. *
117 !* J. Remote Sens., 11, 659-668, 1990. *
118 !* Tanre, D., C. Deroo, P. Duhaut, M. Herman, J. J. Morcrette, J. Perbos,*
119 !* and P. Y. Deschamps, Simulation of the satellite signal in the *
120 !* solar spectrum (5S), Users' Guide (U. S. T. De Lille, 59655 *
121 !* Villeneu d'Ascq, France: Laboratoire d'Optique Atmospherique), *
122 !* 1986. *
123 !* Thuillier, G., et al., Solar irradiance reference spectra for two *
124 !* solar active levels, Adv. Space Res., 34, 256-261, 2004. *
125 !* Vane, G., R. O. Green, T. G. Chrien, H. T. Enmark, E. G. Hansen, and *
126 !* W. M. Porter, The Airborne Visible/Infrared Imaging Spectrometer, *
127 !* Remote Sens. Env., 44, 127-143, 1993. *
128 !* Vane, G. (Ed), Airborne visible/infrared imaging spectrometer *
129 !* (AVIRIS), JPL Publ. 87-38, Jet Propul. Lab, Pasadena, Calif., 1987.*
130 !* Vermote, E., D. Tanre, J. L. Deuze, M. Herman, and J. J. Morcrette, *
131 !* Second simulation of the satellite signal in the solar spectrum *
132 !* (6S), 6S User's Guide Version 1, NASA-GSFC, Greenbelt, Maryland, *
133 !* 134 pages, 1994. *
134 !* *
135 !********************************************************************************
136  ********************************************************************************
137  */
138 
139 #include "atrem_corl1v3.h"
140 #include <sensorInfo.h>
141 #include <timeutils.h>
142 #include <math.h>
143 #define BUFSZ 100
144 #define MINDEGCHANGE 55 // Minimum change (degrees squared) in sum of zenith and azimuth angle squared
145 // for both solar and sensor angles before recalculating tran_table (transmittance table)
146 
147 int get_atrem_cor(l1str *l1rec, int32_t ip, float *rhot, float *tg_tot, float *tg_sol, float *tg_sen) {
148  int32_t sensorID = l1rec->l1file->sensorID;
149 
150  static paramstr P;
151  int i, j, nb, modnum = input->atrem_model;
152  static int firstCall = 1;
153  int32_t nbands = l1rec->l1file->nbands;
154  static double prev_ddeg = -999, prev_max_senz, prev_min_senz;
155  static float **angle_limit, *ang_senz, *ang_solz;
156  static int n_senz, n_solz, so_ang, se_ang;
157  static int prev_modnum;
158  float limitang, *anglelimit;
159 
160  //Initialize input paramters
161 
162 
163  if (firstCall == 1) { //|| prev_dist > MINDEGCHANGE) {
164  //INIT
165  init_atrem(sensorID, &P, l1rec, nbands);
166  prev_modnum = P.model;
167 
168  if (P.dogeom == 0) {
169  printf("ATREM: Reading geometry limiting angles\n");
170  if (get_angle_limits(&anglelimit, &ang_senz, &ang_solz, &n_senz, &n_solz)) {
171  printf("-E- %s line %d : Error reading angle_limit file.\n",
172  __FILE__, __LINE__);
173  exit(FATAL_ERROR);
174  }
175 
176  if ((angle_limit = (float **) calloc(n_senz, sizeof (float *))) == NULL) {
177  printf("-E- : Error allocating memory to tindx\n");
178  exit(FATAL_ERROR);
179  }
180 
181  for (j = 0; j < n_senz; j++) {
182  if ((angle_limit[j] = (float *) calloc(n_solz, sizeof (float))) == NULL) {
183  printf("-E- : Error allocating memory to tindx\n");
184  exit(FATAL_ERROR);
185  }
186  for (i = 0; i < n_solz; i++) {
187  angle_limit[j][i] = *(anglelimit + j * (n_solz) + i);
188  //printf("RJH:2: angle_limit: %f %f %f\n",ang_solz[i],ang_senz[j],angle_limit[j][i]);
189 
190  }
191  }
192 
193  prev_max_senz = l1rec->senz[ip];
194  prev_min_senz = l1rec->senz[ip];
195  }
196  }
197 
198  if (P.dogeom == 0) {
199  if (l1rec->senz[ip] > prev_max_senz) prev_max_senz = l1rec->senz[ip];
200  if (l1rec->senz[ip] < prev_min_senz) prev_min_senz = l1rec->senz[ip];
201 
202  prev_ddeg = fabs(prev_max_senz - prev_min_senz);
203 
204  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);
205  //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]);
206  }
207  if (firstCall == 1 || prev_ddeg >= limitang || P.dogeom != 0) {
208 
209  //printf("Calculating Transmittance table for Atrem correction, ip=%d\n",ip);
210  //start_time = now();
211  //printf("\nBegin get_atrem_cor processing at %s\n\n", ydhmsf(start_time,'L'));
212 
213 
214  geometry_l2gen_.splitpaths = 0;
215  geometry_l2gen_.senzn_l2 = l1rec->senz[ip];
216  geometry_l2gen_.senaz_l2 = l1rec->sena[ip];
217  geometry_l2gen_.solzn_l2 = l1rec->solz[ip];
218  getinput3_.vrto3 = l1rec->oz[ip];
219 
220  //geometry_(); //FORTRAN
221  geometry();
222  //printf("Processed geometry after %6.0f seconds\n",now()-start_time);
223 
224  //init_speccal_(); //FORTRAN
225  init_spectral_calculations();
226  //printf("Processed init_speccal after %6.0f seconds\n",now()-start_time);
227 
228  //tran_table_(); //FORTRAN
229  tran_table();
230  //printf("Processed tran_table after %6.0f seconds\n",now()-start_time);
231 
232  P.nh2o = init_speccal3_.nh2o;
233  P.start_ndx[0] = init_speccal6_.ist1 - 1; // Fortran array index starts at 1, C at 0
234  P.start_ndx[1] = init_speccal6_.ist2 - 1;
235  P.end_ndx[0] = init_speccal6_.ied1 - 1;
236  P.end_ndx[1] = init_speccal6_.ied2 - 1;
237  P.start_p94 = init_speccal6_.istp94 - 1;
238  P.end_p94 = init_speccal6_.iedp94 - 1;
239  P.start_ndx[2] = init_speccal7_.ist3 - 1;
240  P.start_ndx[3] = init_speccal7_.ist4 - 1;
241  P.end_ndx[2] = init_speccal7_.ied3 - 1;
242  P.end_ndx[3] = init_speccal7_.ied4 - 1;
243  P.start_1p14 = init_speccal7_.ist1p14 - 1;
244  P.end_1p14 = init_speccal7_.ied1p14 - 1;
245  P.wt1 = init_speccal8_.wt1;
246  P.wt2 = init_speccal8_.wt2;
247  P.wt3 = init_speccal8_.wt3;
248  P.wt4 = init_speccal8_.wt4;
249  P.start2 = init_speccal10_.istrt2 - 1;
250  P.end2 = init_speccal10_.iend2 - 1;
251  P.ncv2 = init_speccal10_.ncv2;
252  P.finst2 = init_speccal10_.finst2;
253  P.natot = init_speccal11_.natot;
254  P.nbtot = init_speccal11_.nbtot;
255  P.nctot = init_speccal11_.nctot;
256  P.ndtot = init_speccal11_.ndtot;
257  P.g_vap_equiv = geometry3_.g_vap_equiv;
258  P.r0p94 = tran_table1_.r0p94;
259  P.r1p14 = tran_table1_.r1p14;
260  P.vaptot = tran_table1_.vaptot;
261  P.trntbl = tran_table1_.trntbl;
262 
263 
264  // printf("PARAM: nb: %d %d %d %d\n",P.nb1,P.nb2,P.nb3,P.nb4);
265  // printf("PARAM:nb: %d %d\n",P.nbp94,P.nb1p14);
266  // printf("PARAM:nh2o: %d\n",P.nh2o);
267  // printf("PARAM:nbands: %d\n",P.nbands);
268  // printf("PARAM:startndx: %d %d %d %d\n",P.start_ndx[0],P.start_ndx[1],P.start_ndx[2],P.start_ndx[3]);
269  // printf("PARAM:endndx: %d %d %d %d\n",P.end_ndx[0],P.end_ndx[1],P.end_ndx[2],P.end_ndx[3]);
270  // printf("PARAM:start_p94: %d\n",P.start_p94);
271  // printf("PARAM:end_p94: %d\n",P.end_p94);
272  // printf("PARAM:start_1p14: %d\n",P.start_1p14);
273  // printf("PARAM:end_1p14: %d\n",P.end_1p14);
274  // printf("PARAM:%d\n",P.start2);
275  // printf("PARAM:%d\n",P.end2);
276  //
277  // printf("PARAM:ncv2: %d\n",P.ncv2);
278  // printf("PARAM:ntot: %d %d %d %d\n",P.natot,P.nbtot,P.nctot,P.ndtot);
279  // printf("PARAM:wt: %f %f %f %f\n",P.wt1,P.wt2,P.wt3,P.wt4);
280  // printf("PARAM:delta: %f %f\n",P.delta,P.delta2);
281  // printf("PARAM:g_vap_equiv: %f \n",P.g_vap_equiv);
282  // printf("PARAM:ip=%d\n",ip);
283  // printf("PARAM:prev_dist=%f\n",prev_ddeg);
284  firstCall = 0;
285 
286  prev_max_senz = l1rec->senz[ip];
287  prev_min_senz = l1rec->senz[ip];
288  }
289 
290  // if (l1rec->iscan != prevscan) {
291  // prevscan = l1rec->iscan;
292  // }
293  if (modnum == 0) {
294  int16_t year, day;
295  double secs;
296  unix2yds(l1rec->scantime, &year, &day, &secs);
297  P.model = getModelNum(*l1rec->lat, day);
298  if (P.model != prev_modnum) {
299  nb = init_tpvmr_nc(P.model);
300  if (nb <= 0) {
301  printf("-E- %s line %d : Atmospheric data could not be initialized.\n",
302  __FILE__, __LINE__);
303  exit(FATAL_ERROR);
304  }
305  //model_adj_(); //FORTRAN
306  model_adjust();
307  prev_modnum = P.model;
308  }
309  }
310  l1rec->wv[ip] = get_atrem(tg_tot, rhot, &P); //returns water vapor(cm)
311 
312  // Calculate the separate path transmittances (solar and sensor) if selected
313  if (input->atrem_splitpaths) {
314  geometry_l2gen_.splitpaths = 1;
315  geometry_l2gen_.water_vapor = l1rec->wv[ip]; //Misunderstanding. This isn't used.
316  geometry_l2gen_.ja = P.ja;
317  geometry_l2gen_.jb = P.jb;
318  geometry_l2gen_.f1a = P.f1a;
319  geometry_l2gen_.f1b = P.f1b;
320  geometry_l2gen_.f2a = P.f2a;
321  geometry_l2gen_.f2b = P.f2b;
322 
323  // Call geometry with the indices from the tran_table
324  //geometry_(); //FORTRAN
325  geometry();
326 
327  // Don't need to call init_speccal since other trace gas transmittances don't depend on water vapor
328 
329  // recalculate transmittances based on separate solar/sensor paths using the tran_table indices from the first call
330  //tran_table_(); //FORTRAN
331  tran_table();
332 
333  for (i = 0; i < P.nbands; i++) {
334  tg_sol[i] = tran_table_l2gen_.tg_sol[i];
335  tg_sen[i] = tran_table_l2gen_.tg_sen[i];
336  }
337 
338  }
339  //printf("Water vapor[%d]=%f\n",ip,l1rec->wv[ip]);
340  //printf("Processed get_atrem after %6.0f seconds\n",now()-start_time);
341 
342 
343  return (0);
344 }
345 
346 float get_atrem(float *tg_tot, float *rhot, paramstr *P) {
347 
348  double const1 = 0, const2 = 0, const3 = 0, const4 = 0, const5 = 0, const6 = 0;
349  double ratio_94c, ratio_94co, ratio_114c, ratio_114co;
350  int32_t i, ja, jb;
351  float clmwvp;
352  /*
353  Arrays related to look-up table:
354  VAPTOT: TOTAL SUN-SURFACE-SENSOR PATH WATER VAPOR IN UNITS OF CM
355  R0P94 : Tc(0.94 um)/(WT1*Tc(0.86)+WT2*Tc(1.02))
356  R1P14 : Tc(1.14 um)/(WT3*Tc(1.05)+WT4*Tc(1.23))
357 
358  Calculating 3-channel ratios from an observed spectrum, using a
359  look-up table procedure to derive the column amount of water vapor
360  and to find the associated simulated transmittance spectrum.
361  */
362  for (i = P->start_ndx[0]; i <= P->end_ndx[0]; i++) {
363  const1 += rhot[i];
364  }
365  const1 /= P->nb1;
366 
367  for (i = P->start_ndx[1]; i <= P->end_ndx[1]; i++) {
368  const2 += rhot[i];
369  }
370  const2 /= P->nb2;
371 
372  // printf("const2=%f nb2=%d istr2=%d ind2=%d\n",const2,P->nb2,P->start_ndx[1],P->end_ndx[1]);
373 
374  for (i = P->start_p94; i <= P->end_p94; i++) {
375  const3 += rhot[i];
376  }
377  const3 /= P->nbp94;
378 
379  ratio_94co = const3 / ((P->wt1 * const1) + (P->wt2 * const2));
380  ratio_94c = ratio_94co;
381 
382  if (ratio_94co > 1.0) {
383  const1 = 0.0;
384 
385  for (i = P->start_ndx[0]; i <= P->end_ndx[0]; i++) {
386  const1 += (1.0 / rhot[i]);
387  }
388  const1 /= P->nb1;
389 
390  const2 = 0.0;
391  for (i = P->start_ndx[1]; i <= P->end_ndx[1]; i++) {
392  const2 += (1.0 / rhot[i]);
393  }
394  const2 /= P->nb2;
395  const3 = 0.0;
396  for (i = P->start_p94; i <= P->end_p94; i++) {
397  const3 += (1.0 / rhot[i]);
398  }
399  const3 /= P->nbp94;
400 
401  ratio_94c = const3 / ((P->wt1 * const1) + (P->wt2 * const2));
402  }
403 
404  debug_atrem.rp94 = ratio_94c;
405 
406  const4 = 0.0;
407  for (i = P->start_ndx[2]; i <= P->end_ndx[2]; i++) {
408  const4 += rhot[i];
409  }
410  const4 /= P->nb3;
411 
412  const5 = 0.0;
413  for (i = P->start_ndx[3]; i <= P->end_ndx[3]; i++) {
414  const5 += rhot[i];
415  }
416  const5 /= P->nb4;
417 
418  const6 = 0.0;
419  for (i = P->start_1p14; i <= P->end_1p14; i++) {
420  const6 += rhot[i];
421  }
422  const6 /= P->nb1p14;
423 
424  /* DEBUG
425  *
426  */
427  debug_atrem.cst1 = const1;
428  debug_atrem.cst2 = const2;
429  debug_atrem.cst3 = const3;
430  debug_atrem.cst4 = const4;
431  debug_atrem.cst5 = const5;
432  debug_atrem.cst6 = const6;
433 
434  ratio_114co = const6 / ((P->wt3 * const4) + (P->wt4 * const5));
435  ratio_114c = ratio_114co;
436 
437  if (ratio_114co > 1.0) {
438 
439  const4 = 0.0;
440  for (i = P->start_ndx[2]; i <= P->end_ndx[2]; i++) {
441  const4 += (1.0 / rhot[i]);
442  }
443  const4 /= P->nb3;
444  for (i = P->start_ndx[3]; i <= P->end_ndx[3]; i++) {
445  const5 += (1.0 / rhot[i]);
446  }
447  const5 /= P->nb4;
448  const6 = 0.0;
449  for (i = P->start_1p14; i <= P->end_1p14; i++) {
450  const6 += (1.0 / rhot[i]);
451  }
452  const6 /= P->nb1p14;
453  ratio_114c = const6 / ((P->wt3 * const4) + (P->wt4 * const5));
454  }
455 
456  debug_atrem.r1p14 = ratio_114c;
457 
458  double delta, deltab, fja, fjap1, fjb, fjbp1, vaptta, vapttb, specav, spec450;
459  static double *speca, *specb;
460 
461  if (!speca) speca = (double*) allocateMemory(P->nbands * sizeof (double), "speca");
462  if (!specb) specb = (double*) allocateMemory(P->nbands * sizeof (double), "specb");
463 
464  ja = P->nh2o / 2;
465  ja = hunt(P->r0p94, P->nh2o, ratio_94c, ja);
466  // printf("xx[%d]=%f xx[%d]=%f JA=%d\n",ja-1,P->r0p94[ja-1],ja,P->r0p94[ja],ja);
467  if (ja >= 0 && ja < P->nh2o) {
468  delta = P->r0p94[ja + 1] - P->r0p94[ja];
469  fja = (P->r0p94[ja + 1] - ratio_94c) / delta;
470  fjap1 = (ratio_94c - P->r0p94[ja]) / delta;
471  vaptta = fja * P->vaptot[ja] + fjap1 * P->vaptot[ja + 1];
472  if (ratio_94co > 1.0) vaptta = -vaptta;
473  } else {
474  if (ja < 0) vaptta = P->vaptot[ja + 1];
475  if (ja > P->nh2o) vaptta = P->vaptot[ja];
476  }
477 
478  if (ratio_94co <= 1.0) {
479  for (i = 0; i < P->nbands; i++) {
480  if (ja >= 0 && ja < P->nh2o - 1) {
481  speca[i] = fja * P->trntbl[ja][i] + fjap1 * P->trntbl[ja + 1][i];
482  } else {
483  if (ja < 0) speca[i] = P->trntbl[ja + 1][i];
484  if (ja >= P->nh2o - 1) speca[i] = P->trntbl[ja][i];
485  }
486  }
487  }
488 
489  if (ratio_94co > 1.0) {
490  for (i = 0; i < P->nbands; i++) {
491  if (ja >= 0 && ja < P->nh2o - 1) {
492  speca[i] = 1.0 / (fja * P->trntbl[ja][i] + fjap1 * P->trntbl[ja + 1][i]);
493  } else {
494  if (ja < 0) speca[i] = 1.0 / P->trntbl[ja + 1][i];
495  if (ja >= NH2OMAXM1) speca[i] = 1.0 / P->trntbl[ja][i];
496  }
497  }
498  }
499 
500  jb = ja;
501 
502  jb = hunt(&P->r1p14[0], P->nh2o, ratio_114c, jb);
503  //printf("RJH: 1p14 JB=%d\n",jb);
504 
505  debug_atrem.jac = ja;
506  debug_atrem.jbc = jb;
507 
508  if (jb >= 0 && jb < P->nh2o - 1) {
509  deltab = P->r1p14[jb + 1] - P->r1p14[jb];
510  fjb = (P->r1p14[jb + 1] - ratio_114c) / deltab;
511  fjbp1 = (ratio_114c - P->r1p14[jb]) / deltab;
512  vapttb = fjb * P->vaptot[jb] + fjbp1 * P->vaptot[jb + 1];
513  if (ratio_114co > 1.0) vapttb = -vapttb;
514  } else {
515  if (jb < 0) vapttb = P->vaptot[jb + 1];
516  if (jb <= P->nh2o - 1) vapttb = P->vaptot[jb];
517  }
518 
519  if (ratio_114co <= 1.0) {
520  for (i = 0; i < P->nbands; i++) {
521  if (jb >= 0 && jb < P->nh2o - 1) {
522  specb[i] = fjb * P->trntbl[jb][i] + fjbp1 * P->trntbl[jb + 1][i];
523  } else {
524  if (jb < 0) specb[i] = P->trntbl[jb + 1][i];
525  if (jb >= P->nh2o - 1) specb[i] = P->trntbl[jb][i];
526  }
527  }
528  }
529  if (ratio_114co > 1.0) {
530  for (i = 0; i < P->nbands; i++) {
531  if (jb >= 0 && jb < P->nh2o - 1) {
532  specb[i] = 1.0 / (fjb * P->trntbl[jb][i] + fjbp1 * P->trntbl[jb + 1][i]);
533  } else {
534  if (jb < 0) specb[i] = 1.0 / P->trntbl[jb + 1][i];
535  if (jb >= P->nh2o - 1) specb[i] = 1.0 / P->trntbl[jb][i];
536  }
537  }
538  }
539 
540  clmwvp = 0.5 * (vaptta + vapttb) / P->g_vap_equiv;
541  spec450 = 1.0 - 0.5 * (speca[P->idx450] + specb[P->idx450]); // should be 0.0
542 
543  // Derivation of surface reflectances
544 
545  for (i = 0; i < P->nbands; i++) {
546  specav = 0.5 * (speca[i] + specb[i]);
547  tg_tot[i] = specav + spec450;
548  //printf("RJH: Atrem: %d %f YY=%f %f %f\n",i,tg_tot[i],rhot[i],rhot[i]/specav,rhot[i]-rhot[i]/specav);
549  }
550 
551  P->ja = ja;
552  P->jb = jb;
553  P->f1a = fja;
554  P->f2a = fjap1;
555  P->f1b = fjb;
556  P->f2b = fjbp1;
557 
558 
559  /*
560  // Smooth the derived surface reflectance spectra
561  // rhot = yy from atrem fortran code
562  // tg_tot = vc from atrem fortran code
563  int ii,j;
564  double truncv;
565 
566  if (P->delta2 > P->delta) {
567 
568  // * First, replace radiances <= 0 near the spectral overlapping parts of the
569  // * four AVIRIS spectrometers by radiances in the nearby AVIRIS' channels.
570 
571  for (i=P->natot-3; i< P->natot+2; i++) {
572  if (rhot[i] <= 0.0) rhot[i] = rhot[i-1];
573  }
574  for (i=P->nbtot-3; i< P->nbtot+2; i++) {
575  if (rhot[i] <= 0.0) rhot[i] = rhot[i-1];
576  }
577  for (i=P->nctot-3; i< P->nctot+2; i++) {
578  if (rhot[i] <= 0.0) rhot[i] = rhot[i-1];
579  }
580  for (i=P->ndtot-3; i< P->ndtot+2; i++) {
581  if (rhot[i] <= 0.0) rhot[i] = rhot[i-1];
582  }
583  for (i=P->start2-1; i<P->end2; i++){
584  truncv = 0.0;
585  ii = i - P->ncv2 - 1;
586  for (j=ii; j<i+P->ncv2; j++) {
587  truncv += rhot[j]*P->finst2[j-ii+1];
588  }
589  rhot[i] = truncv;
590  }
591  }
592  */
593 
594  return (clmwvp);
595 
596 }
597 
598 /*********************************************************************************
599  * *
600  * Name: HUNT *
601  * Purpose: finds the element in array XX that is closest to value X. Array AA *
602  * must be monotonic, either increasing or decreasing. *
603  * Parameters: XX - array to search *
604  * N - number of elements in the array *
605  * X - element to search for closest match *
606  * JLO - index of the closest matching element *
607  * Algorithm: this subroutine was copied from Numerical Recipes *
608  * Modified for C and cleaned up by *
609  * Richard Healy SAIC/NASA-GSFC 2/19/2015 *
610  * (Who didn't realize it was a Numerical Recipe's function until finished) *
611  * Globals used: none *
612  * Global output: none *
613  * Return Codes: none *
614  * Special Considerations: none *
615  * *
616  **********************************************************************************
617  *
618  */
619 int32_t hunt(float *xx, int32_t n, double x, int32_t jlo) {
620  int32_t inc, jhi, jm;
621  int ascnd;
622 
623  ascnd = (xx[n - 1] >= xx[0]);
624 
625  inc = 1;
626  if (jlo >= 0 && jlo < n) {
627 
628  if ((x >= xx[jlo]) == ascnd) {
629  jhi = jlo + inc;
630  while (jhi < n && ((x >= xx[jhi]) == ascnd)) {
631  jlo = jhi;
632  inc += inc;
633  jhi = jlo + inc;
634  }
635  if (jhi >= n)
636  jhi = n;
637 
638  } else {
639  jhi = jlo;
640  jlo = jhi - inc;
641  while (jlo >= 0 && ((x < xx[jlo]) == ascnd)) {
642  jhi = jlo;
643  inc += inc;
644  jlo = jhi - inc;
645  }
646  if (jlo < 0)
647  jlo = -1;
648 
649  }
650 
651  } else {
652  jlo = 0;
653  jhi = n + 1;
654  }
655  while (jhi - jlo != 1) {
656 
657  jm = (jhi + jlo) / 2;
658 
659  if ((x >= xx[jm]) == ascnd) {
660  jlo = jm;
661  } else {
662  jhi = jm;
663  }
664  }
665 
666  if (x == xx[n - 1]) jlo = n - 2;
667  if (x == xx[0]) jlo = 0;
668 
669  return (jlo);
670 }
671 
672 int init_tpvmr(int model) {
673 
674  int i, nb;
675  printf("ATREM: Initializing ATM for model number = %d\n", model);
676  model--;
677  if (model < 0 || model > MODELMAX) {
678  printf("-E- %sline %d: Invalid ATM Model number\n Value must be between 1 and 7\n: get_atrem_cor3\n",
679  __FILE__, __LINE__);
680  exit(FATAL_ERROR);
681  }
682  nb = tpvmr_init1_.tpvmr[0][model];
683 
684  for (i = 0; i < nb; i++) {
685  getinput3_.h[i] = tpvmr_init1_.tpvmr[1 + (4 * i)][model];
686  //Convert the atmospheric pressure from "mb" to "atm.":
687  getinput3_.p[i] = tpvmr_init1_.tpvmr[2 + (4 * i)][model] / 1013.;
688  getinput3_.t[i] = tpvmr_init1_.tpvmr[3 + (4 * i)][model];
689  //Convert the VMR from the ppm unit in the model to absolute unit
690  getinput3_.vmr[i] = tpvmr_init1_.tpvmr[4 + (4 * i)][model]*1.0E-06;
691  }
692 
693  for (i = nb; i < MODELMAX; i++) {
694  getinput3_.h[i] = 1000.;
695  getinput3_.p[i] = 0.0;
696  getinput3_.t[i] = 300;
697  getinput3_.vmr[i] = 0.0;
698  }
699  getinput3_.nb = nb;
700  getinput3_.nl = nb - 1;
701 
702 
703 
704  return nb;
705 }
706 
707 void get_tpvmr(size_t layers, size_t models, int sds_id,
708  char filename[FILENAME_MAX], char *varname, float* var_a) {
709  size_t start[3], count[3];
710  int retval;
711  int xid;
712  start[0] = 0;
713  start[1] = 0;
714  start[2] = 0;
715  count[0] = models;
716  count[1] = layers;
717  count[2] = 0;
718  retval = nc_inq_varid(sds_id, varname, &xid);
719  if (retval != NC_NOERR) {
720  fprintf(stderr,
721  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
722  __FILE__, __LINE__, filename, varname);
723  exit(FATAL_ERROR);
724  }
725  retval = nc_get_vara_float(sds_id, xid, start, count, var_a);
726  if (retval != NC_NOERR) {
727  fprintf(stderr,
728  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
729  __FILE__, __LINE__, filename, varname);
730  exit(FATAL_ERROR);
731  }
732 }
733 
734 int init_tpvmr_nc(int model) {
735 
736  int i;
737  char *filedir;
738  char filename[FILENAME_MAX];
739  char *tpvmr_file = "atrem_tpvmr.nc";
740  int xid, yid, retval, sds_id;
741  static size_t models, layers;
742  static float *h, *p, *t, *vmr;
743 
744  if (!h) {
745 
746  if ((filedir = getenv("OCDATAROOT")) == NULL) {
747  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
748  return (-1);
749  }
750 
751  strcpy(filename, filedir);
752  strcat(filename, "/common/");
753  strcat(filename, tpvmr_file);
754  strcat(filename, "\0");
755  printf("ATREM: Initializing arrays for Atrem atmospheric models\n");
756 
757  retval = nc_open(filename, NC_NOWRITE, &sds_id);
758  if (retval != NC_NOERR) {
759  fprintf(stderr, "-E- %s line %d: nc_open(%s) failed.\n",
760  __FILE__, __LINE__, filename);
761  exit(FATAL_ERROR);
762  }
763  // Now read in the absorption coefficients
764  retval = nc_inq_dimid(sds_id, "n_layers", &xid) \
765  + nc_inq_dimlen(sds_id, xid, &layers) \
766  + nc_inq_dimid(sds_id, "n_models", &yid) \
767  + nc_inq_dimlen(sds_id, yid, &models);
768 
769  if (retval) {
770  fprintf(stderr, "-E- %s line %d: nc_inq_dimid(%s) failed.%d %d %d \n",
771  __FILE__, __LINE__, filename, (int) layers, (int) models, retval);
772  exit(FATAL_ERROR);
773  }
774 
775  h = (float *) calloc(layers*models, sizeof (float));
776  if (!(h)) {
777  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for h.\n",
778  __FILE__, __LINE__);
779  exit(FATAL_ERROR);
780  }
781  p = (float *) calloc(layers*models, sizeof (float));
782  if (!(p)) {
783  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for p.\n",
784  __FILE__, __LINE__);
785  exit(FATAL_ERROR);
786  }
787  t = (float *) calloc(layers*models, sizeof (float));
788  if (!(t)) {
789  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for t.\n",
790  __FILE__, __LINE__);
791  exit(FATAL_ERROR);
792  }
793  vmr = (float *) calloc(layers*models, sizeof (float));
794  if (!(vmr)) {
795  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for vmr.\n",
796  __FILE__, __LINE__);
797  exit(FATAL_ERROR);
798  }
799  get_tpvmr(layers, models, sds_id, filename, "h", h);
800  get_tpvmr(layers, models, sds_id, filename, "p", p);
801  get_tpvmr(layers, models, sds_id, filename, "t", t);
802  get_tpvmr(layers, models, sds_id, filename, "vmr", vmr);
803  }
804 
805  printf("ATREM: Initializing ATM for model number = %d\n", model);
806  model--;
807 
808  if (model < 0 || model > MODELMAX) {
809  printf("-E- %sline %d: Invalid ATM Model number\n Value must be between 1 and 7\n: get_atrem_cor3\n",
810  __FILE__, __LINE__);
811  exit(FATAL_ERROR);
812  }
813 
814  for (i = 0; i < layers; i++) {
815  getinput3_.h[i] = h[model * layers + i];
816  //Convert the atmospheric pressure from "mb" to "atm.":
817  getinput3_.p[i] = p[model * layers + i] / 1013.;
818  getinput3_.t[i] = t[model * layers + i] / 1013.;
819  //Convert the VMR from the ppm unit in the model to absolute unit
820  getinput3_.vmr[i] = vmr[model * layers + i]*1.0e-06;
821  }
822 
823  for (i = layers; i < MODELMAX; i++) {
824  getinput3_.h[i] = 1000.;
825  getinput3_.p[i] = 0.0;
826  getinput3_.t[i] = 300;
827  getinput3_.vmr[i] = 0.0;
828  }
829  getinput3_.nb = layers;
830  getinput3_.nl = layers - 1;
831 
832 
833 
834  return layers;
835 }
836 
837 int getModelNum(float lat, int32_t day) {
838  // Determine which atmospheric model to use
839  // 1 = tropical
840  // 2 = mid latitude summer
841  // 3 = mid latitude winter
842  // 4 = subarctic summer
843  // 5 = subarctic winter
844  // 6 = US standard 1962
845  int16 mon = (int) day / 31 + 1; // month of year (no need for perfection..at least according to the sea surface salinity reference algorithm)
846 
847  if (fabs(lat) < 30) return (1);
848  else {
849  switch (mon) {
850  case 12:
851  case 1:
852  case 2:
853  if (lat < 60 && lat > 30) return (3);
854  if (lat < -30 && lat >-60) return (2);
855  if (lat > 60) return (5);
856  return (4);
857  break;
858  case 6:
859  case 7:
860  case 8:
861  if (lat < 60 && lat > 30) return (2);
862  if (lat < -30 && lat >-60) return (3);
863  if (lat > 60) return (4);
864  return (5);
865  break;
866  default:
867  return (6);
868 
869  }
870 
871 
872  }
873 
874 
875 }
876 
877 int init_atrem(int32_t sensorID, paramstr *P, l1str *l1rec, int32_t nbands) {
878 
879  int32_t nb, atrem_opt = input->atrem_opt, atrem_full = input->atrem_full, atrem_geom = input->atrem_geom;
880  int32_t atrem_splitpaths = input->atrem_splitpaths, atrem_model = input->atrem_model, gas_opt = input->gas_opt;
881  float *fwhm, *xppp;
882  float nwave;
883  int i;
884  char *filedir;
885  char filename[FILENAME_MAX];
886  int *model, flag_gas;
887 
888  model = (int *) calloc(1, sizeof (int));
889  fwhm = (float *) calloc(1, sizeof (float));
890  xppp = (float *) calloc(1, sizeof (float));
891 
892  if ((filedir = getenv("OCDATAROOT")) == NULL) {
893  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
894  return (-1);
895  }
896 
897  strcpy(filename, filedir);
898  strcat(filename, "/");
899  strcat(filename, sensorId2SensorDir(sensorID));
900  strcat(filename, "/\0");
901 
902  //INIT - get_input_ must be called before any common block structures are filled in.
903  // It will also fill some of the parameters with default values
904  //get_input_();
905 
906  strcpy(input_l2gen_.filename, filename);
907  input_l2gen_.dln = strlen(filename);
908  //printf("filedir=%s Filename=%s len=%d\n",filedir, filename,input_l2gen_.dln);
909 
910  if (!(fwhm = l1rec->l1file->fwhm)) {
911  nwave = rdatreminfo(sensorID, 0, "fwhm", (void **) &fwhm);
912  }
913 
914  getinput4_.fwhm = (float*) allocateMemory(nbands * sizeof (float), "fwhm");
915  getinput4_.wavobs = (float*) allocateMemory(nbands * sizeof (float), "wavobs");
916 
917  for (i = 0; i < nbands; i++) {
918  getinput4_.wavobs[i] = l1rec->l1file->fwave[i] / 1000.; // Convert nm to microns
919  getinput4_.fwhm[i] = fwhm[i]; //fwhm should already be in microns
920  //printf("->RJH:fwave(%d)=%f fwhm(%d) = %f \n",i+1,getinput4_.wavobs[i],i,fwhm[i]);
921  }
922 
923  P->idx450 = bindex_get(450); // get the index nearest the 450nm wavelength
924  if (P->idx450 < 0) {
925  printf("-E- %s line %d : Unable to determine 450 nm index from spectrum.\n",
926  __FILE__, __LINE__);
927  exit(FATAL_ERROR);
928  }
929  // for (i=0;i<nbands;i++) {
930  // getinput4_.fwhm[i] = fwhm[i];
931  // // printf("->RJH:fwave(%d)=%f fwhm(%d) = %f\n",i+1,getinput4_.wavobs[i],i,fwhm[i]);
932  // }
933  /*
934  // Determine which atmospheric model to use
935  // call to get_input sets h,p,t, and vmr arrays to default model 1
936  // 1 = tropical
937  // 2 = mid latitude summer
938  // 3 = mid latitude winter
939  // 4 = subarctic summer
940  // 5 = subarctic winter
941  // 6 = US standard 1962
942  // 7 = user defined model
943 
944  //INIT
945  //Get Atmospheric model to use
946  */
947  *model = 6;
948 
949  if (atrem_model == 0) {
950  int16_t year, day;
951  double secs;
952  unix2yds(l1rec->scantime, &year, &day, &secs);
953  *model = getModelNum(*l1rec->lat, day);
954  } else {
955  if (atrem_model <= 6) {
956  *model = atrem_model;
957  // nwave = rdatreminfo(sensorID, 0, "model", (void **)&model);
958  } else {
959  printf("-E- %s line %d : Invalid atmospheric model, atrem_model = %d. Valid range is 0-6.\n",
960  __FILE__, __LINE__, atrem_model);
961  exit(FATAL_ERROR);
962  }
963  }
964  P->model = *model;
965  nb = init_tpvmr_nc(P->model);
966 
967  if (nb <= 0) {
968  printf("-E- %s line %d : Atmospheric data could not be initialized.\n",
969  __FILE__, __LINE__);
970  exit(FATAL_ERROR);
971  }
972 
973  getinput5_.nbands = nbands; //l1rec->nbands;
974  getinput5_.dlt = 0; //used internally to the fortran routines
975  getinput5_.dlt2 = 0; //to determine what to use for the width of the bands
976  getinput5_.hsurf = 0; //mean surface elevation
977  getinput14_.xppp = 700; //Plane/Satellite altitude (km)
978 
979  if (getinput5_.hsurf < 0 || getinput5_.hsurf > getinput3_.h[nb - 2]) {
980  printf("-E- %s line %d : Elevation=%f must be less than the maximum elevation in the atmospheric model.\n",
981  __FILE__, __LINE__, getinput5_.hsurf);
982  exit(FATAL_ERROR);
983  }
984 
985  nwave = rdatreminfo(sensorID, 0, "xppp", (void **) &xppp);
986 
987  if (getinput14_.xppp < getinput5_.hsurf) {
988  printf("-E- %s line %d : Sensor altitude=%f must be greater than the bottom surface elevation.\n",
989  __FILE__, __LINE__, getinput14_.xppp);
990  exit(FATAL_ERROR);
991 
992  }
993  /*
994  // Get ranges on curve for the two atmospheric windows surrounding the 1.14-um
995  // water vapor absorption feature and the center point of the 1.14-um water
996  // vapor absorption feature. Enter:
997  // 1. the midpoint of third window (0.6-2.5)
998  // 2. number of points to average for third window (1-10)
999  // 3. the midpoint of fourth window (0.6-2.5)
1000  // 4. number of points to average for fourth window (1-10)
1001  // 5. the midpoint of 1.14-um absorption feature (0.6-2.5)
1002  // 6. the number of points to average for the absorption feature (1-30)
1003  */
1004  // This is the default for HICO HS
1005  float *wndow1, *wndow2, *wp94c, *wndow3, *wndow4, *w1p14c;
1006 
1007  wndow1 = (float *) calloc(1, sizeof (float));
1008  wndow2 = (float *) calloc(1, sizeof (float));
1009  wndow3 = (float *) calloc(1, sizeof (float));
1010  wndow4 = (float *) calloc(1, sizeof (float));
1011  wp94c = (float *) calloc(1, sizeof (float));
1012  w1p14c = (float *) calloc(1, sizeof (float));
1013 
1014  getinput6_.wndow1 = 0.705;
1015  getinput6_.wndow2 = 0.745;
1016  getinput6_.wp94c = 0.725;
1017  getinput6_.wndow3 = 0.805;
1018  getinput6_.wndow4 = 0.845;
1019  getinput6_.w1p14c = 0.825;
1020 
1021  nwave = rdatreminfo(sensorID, 0, "window1", (void **) &wndow1);
1022  nwave = rdatreminfo(sensorID, 0, "window2", (void **) &wndow2);
1023  nwave = rdatreminfo(sensorID, 0, "window3", (void **) &wndow3);
1024  nwave = rdatreminfo(sensorID, 0, "window4", (void **) &wndow4);
1025  nwave = rdatreminfo(sensorID, 0, "wp94c", (void **) &wp94c);
1026  nwave = rdatreminfo(sensorID, 0, "w1p14c", (void **) &w1p14c);
1027 
1028  getinput6_.wndow1 = *wndow1;
1029  getinput6_.wndow2 = *wndow2;
1030  getinput6_.wp94c = *wp94c;
1031  getinput6_.wndow3 = *wndow3;
1032  getinput6_.wndow4 = *wndow4;
1033  getinput6_.w1p14c = *w1p14c;
1034 
1035  if (getinput6_.wndow1 < 0.6 || getinput6_.wndow1 > 2.5) {
1036  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
1037  "in the .94-um region1:%f\nValid values are 0.6-2.5\n", getinput6_.wndow1);
1038  printf("-E- %s line %d : See stderr output for details.\n",
1039  __FILE__, __LINE__);
1040  exit(FATAL_ERROR);
1041  }
1042  if (getinput6_.wndow2 < 0.6 || getinput6_.wndow2 > 2.5) {
1043  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
1044  "in the .94-um region2:%f\nValid values are 0.6-2.5\n", getinput6_.wndow2);
1045  printf("-E- %s line %d : See stderr output for details.\n",
1046  __FILE__, __LINE__);
1047  exit(FATAL_ERROR);
1048  }
1049  if (getinput6_.wp94c <= getinput6_.wndow1 || getinput6_.wp94c >= getinput6_.wndow2) {
1050  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
1051  "in the .94-um region:%f\nValid range is: %f < value < %f \n", getinput6_.wp94c, getinput6_.wndow1, getinput6_.wndow2);
1052  printf("-E- %s line %d : See stderr output for details.\n",
1053  __FILE__, __LINE__);
1054  exit(FATAL_ERROR);
1055  }
1056  if (getinput6_.wndow3 < 0.6 || getinput6_.wndow3 > 2.5) {
1057  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
1058  "in the .94-um region3:%f\nValid values are 0.6-2.5\n", getinput6_.wndow3);
1059  printf("-E- %s line %d : See stderr output for details.\n",
1060  __FILE__, __LINE__);
1061  exit(FATAL_ERROR);
1062  }
1063  if (getinput6_.wndow4 < 0.6 || getinput6_.wndow4 > 2.5) {
1064  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
1065  "in the .94-um region4:%f\nValid values are 0.6-2.5\n", getinput6_.wndow4);
1066  printf("-E- %s line %d : See stderr output for details.\n",
1067  __FILE__, __LINE__);
1068  exit(FATAL_ERROR);
1069  }
1070  if (getinput6_.w1p14c <= getinput6_.wndow3 || getinput6_.w1p14c >= getinput6_.wndow4) {
1071  fprintf(stderr, "Invalid wavelength position for first atmospheric window " \
1072  "in the .94-um region:%f\nValid range is: %f < value < %f \n", getinput6_.w1p14c, getinput6_.wndow3, getinput6_.wndow4);
1073  printf("-E- %s line %d : See stderr output for details.\n",
1074  __FILE__, __LINE__);
1075  exit(FATAL_ERROR);
1076  }
1077 
1078  // Set full_calc to 1 to turn on explicit, full calculation of
1079  // water vapor correction. Setting value to 0 turns on speedy
1080  // h2o absorption calculation and sets all other gas options to 0
1081 
1082  // int32_t *full_calc;
1083  // full_calc = (int32_t *) calloc(1,sizeof(int32_t));
1084  // if (atrem_full != 0) {
1085  // *full_calc = (atrem_full < 0? 0:1);
1086  // } else {
1087  // nwave = rdatreminfo(sensorID, 0, "full_calc", (void **) &full_calc);
1088  // }
1089 
1090  /*
1091  * Allocate memory to global variables
1092  */
1093 
1094  init_speccal1_.tran_hi_others = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_others");
1095  init_speccal12_.wavln_med = (float*) allocateMemory(NP_MED * sizeof (float), "wavln_med");
1096  init_speccal12_.wavln_std = (float*) allocateMemory(NP_STD * sizeof (float), "wavln_med");
1097  init_speccal13_.index_med = (int32_t*) allocateMemory(NP_MED * sizeof (int32_t), "index_med");
1098  init_speccal13_.wavln_med_index = (float*) allocateMemory(NP_MED * sizeof (float), "wavln_med_index");
1099  init_speccal15_.ncvhf = (int32_t*) allocateMemory(nbands * sizeof (int32_t), "ncvhf");
1100 
1101  init_speccal13_.tran_med_index = (float **) malloc(NH2OMAX * sizeof (float *));
1102  for (i = 0; i < NH2OMAX; i++)
1103  init_speccal13_.tran_med_index[i] = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index");
1104 
1105  init_speccal15_.finstr = (float **) malloc(nbands * sizeof (float *));
1106  for (i = 0; i < nbands; i++)
1107  init_speccal15_.finstr[i] = (float*) allocateMemory(NINSTR_MAX * sizeof (float), "finstr");
1108 
1109  init_speccal16_.tran_o3_std = (float*) allocateMemory(NO3PT * sizeof (float), "tran_o3_std");
1110  init_speccal17_.tran_no2_std = (float*) allocateMemory(NO3PT * sizeof (float), "tran_no2_std");
1111 
1112  tran_table1_.vaptot = (float*) allocateMemory(NH2OMAX * sizeof (float), "vaptot");
1113  tran_table1_.r0p94 = (float*) allocateMemory(NH2OMAX * sizeof (float), "r0p94");
1114  tran_table1_.r1p14 = (float*) allocateMemory(NH2OMAX * sizeof (float), "r1p14");
1115  tran_table1_.trntblo = (float*) allocateMemory(nbands * sizeof (float), "trntblo");
1116 
1117  tran_table1_.trntbl = (float **) malloc(NH2OMAX * sizeof (float *));
1118  tran_table1_.tran_kd = (float **) malloc(NH2OMAX * sizeof (float *));
1119  tran_table1_.diff_tran = (float **) malloc(NH2OMAX * sizeof (float *));
1120 
1121  for (i = 0; i < NH2OMAX; i++) {
1122  tran_table1_.trntbl[i] = (float*) allocateMemory(nbands * sizeof (float), "trntbl");
1123  tran_table1_.tran_kd[i] = (float*) allocateMemory(nbands * sizeof (float), "tran_kd");
1124  tran_table1_.diff_tran[i] = (float*) allocateMemory(nbands * sizeof (float), "diff_tran");
1125  }
1126 
1127  tran_table_l2gen_.tg_sol = (float*) allocateMemory(nbands * sizeof (float), "tg_sol");
1128  tran_table_l2gen_.tg_sen = (float*) allocateMemory(nbands * sizeof (float), "tg_sen");
1129  tran_table_l2gen_.tg_solo = (float*) allocateMemory(nbands * sizeof (float), "tg_solo");
1130  tran_table_l2gen_.tg_seno = (float*) allocateMemory(nbands * sizeof (float), "tg_seno");
1131 
1132 
1133  if (atrem_splitpaths) {
1134 
1135  //allocate memory
1136 
1137  tran_tables_.tran_hi_sa [0] = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_sa");
1138  tran_tables_.tran_hi_sa [1] = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_sa");
1139  tran_tables_.tran_hi_sap1[0] = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_sap1");
1140  tran_tables_.tran_hi_sap1[1] = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_sap1");
1141  tran_tables_.tran_hi_sb [0] = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_sb");
1142  tran_tables_.tran_hi_sb [1] = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_sb");
1143  tran_tables_.tran_hi_sbp1[0] = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_sbp1");
1144  tran_tables_.tran_hi_sbp1[1] = (float*) allocateMemory(NP_HI * sizeof (float), "tran_hi_sbp1");
1145 
1146  tran_tables1_.tran_med_index_sa_sol = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index_sa_sol");
1147  tran_tables1_.tran_med_index_sa_sen = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index_sa_sen");
1148  tran_tables1_.tran_med_index_sap1_sol = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index_sap1_sol");
1149  tran_tables1_.tran_med_index_sap1_sen = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index_sap1_sen");
1150  tran_tables1_.tran_med_index_sb_sol = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index_sb_sol");
1151  tran_tables1_.tran_med_index_sb_sen = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index_sb_sen");
1152  tran_tables1_.tran_med_index_sbp1_sol = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index_sbp1_sol");
1153  tran_tables1_.tran_med_index_sbp1_sen = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_index_sbp1_sen");
1154 
1155  tran_tables1_.tran_med_sa_sol = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_sa_sol");
1156  tran_tables1_.tran_med_sa_sen = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_sa_sen");
1157  tran_tables1_.tran_med_sap1_sol = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_sap1_sol");
1158  tran_tables1_.tran_med_sap1_sen = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_sap1_sen");
1159  tran_tables1_.tran_med_sb_sol = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_sb_sol");
1160  tran_tables1_.tran_med_sb_sen = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_sb_sen");
1161  tran_tables1_.tran_med_sbp1_sol = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_sbp1_sol");
1162  tran_tables1_.tran_med_sbp1_sen = (float*) allocateMemory(NP_MED * sizeof (float), "tran_med_sbp1_sen");
1163 
1164  tran_tables1_.tran_std_sa_sol = (float*) allocateMemory(NP_STD * sizeof (float), "tran_std_sa_sol");
1165  tran_tables1_.tran_std_sa_sen = (float*) allocateMemory(NP_STD * sizeof (float), "tran_std_sa_sen");
1166  tran_tables1_.tran_std_sap1_sol = (float*) allocateMemory(NP_STD * sizeof (float), "tran_std_sap1_sol");
1167  tran_tables1_.tran_std_sap1_sen = (float*) allocateMemory(NP_STD * sizeof (float), "tran_std_sap1_sen");
1168  tran_tables1_.tran_std_sb_sol = (float*) allocateMemory(NP_STD * sizeof (float), "tran_std_sb_sol");
1169  tran_tables1_.tran_std_sb_sen = (float*) allocateMemory(NP_STD * sizeof (float), "tran_std_sb_sen");
1170  tran_tables1_.tran_std_sbp1_sol = (float*) allocateMemory(NP_STD * sizeof (float), "tran_std_sbp1_sol");
1171  tran_tables1_.tran_std_sbp1_sen = (float*) allocateMemory(NP_STD * sizeof (float), "tran_std_sbp1_sen");
1172 
1173  if (atrem_full == 0) {
1174  printf("ATREM: Turning on atrem_full because you selected atrem_splitpaths\n");
1175  atrem_full = 1;
1176  }
1177  }
1178  if (atrem_full > 0)
1179  printf("ATREM: Warning : full_calc !=0. Atrem will calculate transmittance table for every pixel\n");
1180 
1181  // getinput5_.full_calc = *full_calc;
1182 
1183  getinput5_.full_calc = atrem_full;
1184 
1185  // int32_t *dogeom;
1186  // dogeom = (int32_t *) calloc(1,sizeof(int32_t));
1187  // if (atrem_geom != 0) {
1188  // *dogeom = (atrem_geom < 0? 0:1);
1189  // } else {
1190  // nwave = rdatreminfo(sensorID, 0, "dogeom", (void **) &dogeom);
1191  // }
1192  //
1193  // P->dogeom = *dogeom;
1194 
1195  if (atrem_geom > 0)
1196  printf("ATREM: Warning : dogeom !=0. Geometry will be calculated every pixel\n");
1197 
1198  P->dogeom = atrem_geom;
1199  //Default for HICO HS
1200  int32_t *nb1, *nb2, *nb3, *nb4, *nbp94, *nb1p14;
1201  nb1 = (int32_t *) calloc(1, sizeof (int32_t));
1202  nb2 = (int32_t *) calloc(1, sizeof (int32_t));
1203  nb3 = (int32_t *) calloc(1, sizeof (int32_t));
1204  nb4 = (int32_t *) calloc(1, sizeof (int32_t));
1205  nbp94 = (int32_t *) calloc(1, sizeof (int32_t));
1206  nb1p14 = (int32_t *) calloc(1, sizeof (int32_t));
1207 
1208  getinput7_.nb1 = 3;
1209  getinput7_.nb2 = 3;
1210  getinput7_.nb3 = 3;
1211  getinput7_.nb4 = 3;
1212  getinput7_.nbp94 = 5;
1213  getinput7_.nb1p14 = 5;
1214 
1215  nwave = rdatreminfo(sensorID, 0, "nb1", (void **) &nb1);
1216  nwave = rdatreminfo(sensorID, 0, "nb2", (void **) &nb2);
1217  nwave = rdatreminfo(sensorID, 0, "nb3", (void **) &nb3);
1218  nwave = rdatreminfo(sensorID, 0, "nb4", (void **) &nb4);
1219  nwave = rdatreminfo(sensorID, 0, "nbp94", (void **) &nbp94);
1220  nwave = rdatreminfo(sensorID, 0, "nb1p14", (void **) &nb1p14);
1221 
1222  getinput7_.nb1 = *nb1;
1223  getinput7_.nb2 = *nb2;
1224  getinput7_.nb3 = *nb3;
1225  getinput7_.nb4 = *nb4;
1226  getinput7_.nbp94 = *nbp94;
1227  getinput7_.nb1p14 = *nb1p14;
1228 
1229  if (getinput7_.nb1 < 1 || getinput7_.nb1 > 50) {
1230  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1231  "the .94-um region:%d\nValid values are 1-50\n", getinput7_.nb1);
1232  printf("-E- %s line %d : See stderr output for details.\n",
1233  __FILE__, __LINE__);
1234  exit(FATAL_ERROR);
1235  }
1236 
1237  if (getinput7_.nb2 < 1 || getinput7_.nb2 > 50) {
1238  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1239  "the .94-um region:%d\nValid values are 1-50\n", getinput7_.nb2);
1240  printf("-E- %s line %d : See stderr output for details.\n",
1241  __FILE__, __LINE__);
1242  exit(FATAL_ERROR);
1243  }
1244  if (getinput7_.nbp94 < 1 || getinput7_.nbp94 > 90) {
1245  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1246  "the .94-um region:%d\nValid values are 1-90\n", getinput7_.nbp94);
1247  printf("-E- %s line %d : See stderr output for details.\n",
1248  __FILE__, __LINE__);
1249  exit(FATAL_ERROR);
1250  }
1251  if (getinput7_.nb3 < 1 || getinput7_.nb3 > 50) {
1252  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1253  "the .94-um region:%d\nValid values are 1-50\n", getinput7_.nb3);
1254  printf("-E- %s line %d : See stderr output for details.\n",
1255  __FILE__, __LINE__);
1256  exit(FATAL_ERROR);
1257  }
1258  if (getinput7_.nb4 < 1 || getinput7_.nb4 > 50) {
1259  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1260  "the .94-um region:%d\nValid values are 1-50\n", getinput7_.nb4);
1261  printf("-E- %s line %d : See stderr output for details.\n",
1262  __FILE__, __LINE__);
1263  exit(FATAL_ERROR);
1264  }
1265  if (getinput7_.nb1p14 < 1 || getinput7_.nb1p14 > 110) {
1266  fprintf(stderr, "Invalid number of channels for first wavelength position in " \
1267  "the .94-um region:%d\nValid values are 1-90\n", getinput7_.nbp94);
1268  printf("-E- %s line %d : See stderr output for details.\n",
1269  __FILE__, __LINE__);
1270  exit(FATAL_ERROR);
1271  }
1272 
1273  //INIT
1274  int32_t *h2o, *co2, *o3, *n2o, *co, *ch4, *o2, *no2;
1275  h2o = (int32_t *) calloc(1, sizeof (int32_t));
1276  co2 = (int32_t *) calloc(1, sizeof (int32_t));
1277  o3 = (int32_t *) calloc(1, sizeof (int32_t));
1278  n2o = (int32_t *) calloc(1, sizeof (int32_t));
1279  co = (int32_t *) calloc(1, sizeof (int32_t));
1280  ch4 = (int32_t *) calloc(1, sizeof (int32_t));
1281  o2 = (int32_t *) calloc(1, sizeof (int32_t));
1282  no2 = (int32_t *) calloc(1, sizeof (int32_t));
1283 
1284  *h2o = 1;
1285 
1286  getinput1_.h2o = 1;
1287  getinput1_.co2 = 0;
1288  getinput1_.o3 = 0;
1289  getinput1_.n2o = 0;
1290  getinput1_.co = 0;
1291  getinput1_.ch4 = 0;
1292  getinput1_.o2 = 0;
1293  getinput1_.no2 = 0;
1294 
1295  // nwave = rdatreminfo(sensorID, 0, "h2o", (void **) &h2o);
1296  nwave = rdatreminfo(sensorID, 0, "co2", (void **) &co2);
1297  nwave = rdatreminfo(sensorID, 0, "o3", (void **) &o3);
1298  nwave = rdatreminfo(sensorID, 0, "n2o", (void **) &n2o);
1299  nwave = rdatreminfo(sensorID, 0, "co", (void **) &co);
1300  nwave = rdatreminfo(sensorID, 0, "ch4", (void **) &ch4);
1301  nwave = rdatreminfo(sensorID, 0, "o2", (void **) &o2);
1302  nwave = rdatreminfo(sensorID, 0, "no2", (void **) &no2);
1303 
1304  getinput1_.co2 = *co2;
1305  getinput1_.o3 = *o3;
1306  getinput1_.n2o = *n2o;
1307  getinput1_.co = *co;
1308  getinput1_.ch4 = *ch4;
1309  getinput1_.o2 = *o2;
1310  getinput1_.no2 = *no2;
1311 
1312  // The typical ozone amount is: 0.28-0.55 (atm-cm). The built-in NO2
1313  // column amount is 5.0E+15 molecules/cm^2.
1314  if (atrem_opt > 0) {
1315  getinput1_.co2 = (atrem_opt & ATREM_CO2) > 0;
1316  getinput1_.o3 = (atrem_opt & ATREM_O3) > 0;
1317  getinput1_.n2o = (atrem_opt & ATREM_N2O) > 0;
1318  getinput1_.co = (atrem_opt & ATREM_CO) > 0;
1319  getinput1_.ch4 = (atrem_opt & ATREM_CH4) > 0;
1320  getinput1_.o2 = (atrem_opt & ATREM_O2) > 0;
1321  getinput1_.no2 = (atrem_opt & ATREM_NO2) > 0;
1322 
1323  }
1324 
1325  printf("ATREM: Gas Options:H2O:1 CO2:%d O3:%d N2O:%d CO:%d CH4:%d O2:%d NO2:%d\n",
1326  getinput1_.co2, getinput1_.o3, getinput1_.n2o, getinput1_.co, getinput1_.ch4, getinput1_.o2, getinput1_.no2);
1327 
1328  // Now check to make sure the same gas option isn't selected from the l2gen gas option.
1329 
1330  flag_gas = 0;
1331  if (gas_opt & H2O_BIT) {
1332  printf("ATREM: cannot be used with gas_opt bit mask=%d (H2O)\n", H2O_BIT);
1333  flag_gas = 1;
1334  }
1335  if (getinput1_.no2 == 1 && (gas_opt & NO2_BIT)) {
1336  printf("ATREM: cannot be used with gas_opt bit mask=%d (NO2)\n", NO2_BIT);
1337  flag_gas = 1;
1338 
1339  }
1340  if (getinput1_.co2 == 1 && (gas_opt & CO2_BIT)) {
1341  printf("ATREM: cannot be used with gas_opt bit mask=%d (CO2)\n", CO2_BIT);
1342  flag_gas = 1;
1343  }
1344  if (getinput1_.o3 == 1 && (gas_opt & O3_BIT)) {
1345  printf("ATREM: cannot be used with gas_opt bit mask=%d (O3)\n", O3_BIT);
1346  flag_gas = 1;
1347  }
1348 
1349  if (flag_gas) {
1350  printf("Error: Conflict using ATREM (gas_opt=16 bitmask) with atrem_opt=%d and gas_opt=%d. " \
1351  "\nPlease resolve. Refer to command line options for atrem_opt and gas_opt.\n", atrem_opt, gas_opt);
1352  exit(1);
1353  }
1354 
1355  float *vrto3, *sno2;
1356  vrto3 = (float *) calloc(1, sizeof (float));
1357  sno2 = (float *) calloc(1, sizeof (float));
1358 
1359  getinput3_.vrto3 = 0.34; //total column ozone amount (atm-cm)
1360  getinput3_.sno2 = 1.0; //NO2 scaling factor (to 5.E+15 molecules/cm^2)
1361 
1362  nwave = rdatreminfo(sensorID, 0, "vrto3", (void **) &vrto3);
1363  nwave = rdatreminfo(sensorID, 0, "sno2", (void **) &sno2);
1364 
1365  getinput3_.vrto3 = *vrto3; //total column ozone amount (atm-cm)
1366  getinput3_.sno2 = *sno2; //NO2 scaling factor (to 5.E+15 molecules/cm^2)
1367 
1368  //model_adj_(); //FORTRAN
1369  model_adjust();
1370 
1371  P->nb1 = getinput7_.nb1;
1372  P->nb2 = getinput7_.nb2;
1373  P->nb3 = getinput7_.nb3;
1374  P->nb4 = getinput7_.nb4;
1375  P->nbp94 = getinput7_.nbp94;
1376  P->nb1p14 = getinput7_.nb1p14;
1377  P->nbands = getinput5_.nbands;
1378  P->delta = getinput5_.dlt;
1379  P->delta2 = getinput5_.dlt2;
1380 
1381  return nwave;
1382 }
1383 
1384 int get_angle_limits(float **anglelimit, float **insenz, float **insolz, int *n_senz, int *n_solz) {
1385 
1386  char filename[FILENAME_MAX];
1387  char *infile = "atrem_angle_limit.nc";
1388  char *filedir;
1389 
1390  char name[H4_MAX_NC_NAME];
1391  char sdsname[H4_MAX_NC_NAME];
1392  int ncid;
1393  int32 sds_id;
1394  int status;
1395  nc_type rh_type; /* variable type */
1396  int dimids[H4_MAX_VAR_DIMS]; /* dimension IDs */
1397  int ndims;
1398  int natts; /* number of attributes */
1399  size_t length;
1400 
1401  static float *senz, *solz;
1402  static float *angle_limit;
1403 
1404  if ((filedir = getenv("OCDATAROOT")) == NULL) {
1405  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
1406  return (-1);
1407  }
1408  strcpy(filename, filedir);
1409  strcat(filename, "/common/");
1410  strcat(filename, infile);
1411  strcat(filename, "\0");
1412 
1413  printf("ATREM: Reading Angle_Limit FILE=%s\n", filename);
1414 
1415  if (nc_open(filename, NC_NOWRITE, &ncid) == NC_NOERR) {
1416 
1417  strcpy(sdsname, "senz");
1418 
1419  status = nc_inq_varid(ncid, sdsname, &sds_id);
1420  if (status != NC_NOERR) {
1421  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1422  __FILE__, __LINE__, sdsname, infile);
1423  exit(1);
1424  }
1425 
1426  status = nc_inq_var(ncid, sds_id, 0, &rh_type, &ndims, dimids,
1427  &natts);
1428 
1429  if (nc_inq_dimlen(ncid, dimids[0], &length) != NC_NOERR) {
1430  nc_inq_dim(ncid, dimids[0], name, &length);
1431  fprintf(stderr,
1432  "-E- %s line %d: could not get size of dimension \"%s\" in netCDF File.\n",
1433  __FILE__, __LINE__, name);
1434  exit(1);
1435  }
1436 
1437  *n_senz = length;
1438 
1439  if ((senz = (float *) calloc(*n_senz, sizeof (float))) == NULL) {
1440  printf("-E- : Error allocating memory to senz\n");
1441  exit(FATAL_ERROR);
1442  }
1443 
1444  if (nc_get_var(ncid, sds_id, senz) != NC_NOERR) {
1445  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1446  __FILE__, __LINE__, sdsname, infile);
1447  exit(1);
1448  }
1449 
1450  strcpy(sdsname, "solz");
1451 
1452  status = nc_inq_varid(ncid, sdsname, &sds_id);
1453  if (status != NC_NOERR) {
1454  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1455  __FILE__, __LINE__, sdsname, infile);
1456  exit(1);
1457  }
1458 
1459  status = nc_inq_var(ncid, sds_id, 0, &rh_type, &ndims, dimids,
1460  &natts);
1461  if (status != NC_NOERR) {
1462  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1463  __FILE__, __LINE__, sdsname, infile);
1464  exit(1);
1465  }
1466 
1467  if (nc_inq_dimlen(ncid, dimids[0], &length) != NC_NOERR) {
1468  nc_inq_dim(ncid, dimids[0], name, &length);
1469  fprintf(stderr,
1470  "-E- %s line %d: could not get size of dimension \"%s\" in netCDF File.\n",
1471  __FILE__, __LINE__, name);
1472  exit(1);
1473  }
1474 
1475  *n_solz = length;
1476 
1477  if ((solz = (float *) calloc(*n_solz, sizeof (float))) == NULL) {
1478  printf("-E- : Error allocating memory to solz\n");
1479  exit(FATAL_ERROR);
1480  }
1481 
1482  if (nc_get_var(ncid, sds_id, solz) != NC_NOERR) {
1483  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1484  __FILE__, __LINE__, sdsname, infile);
1485  exit(1);
1486  }
1487 
1488  strcpy(sdsname, "angle_limit");
1489 
1490  status = nc_inq_varid(ncid, sdsname, &sds_id);
1491  if (status != NC_NOERR) {
1492  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1493  __FILE__, __LINE__, sdsname, infile);
1494  exit(1);
1495  }
1496 
1497  status = nc_inq_var(ncid, sds_id, 0, &rh_type, &ndims, dimids,
1498  &natts);
1499  if (status != NC_NOERR) {
1500  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1501  __FILE__, __LINE__, sdsname, infile);
1502  exit(1);
1503  }
1504  if (ndims != 2) {
1505  fprintf(stderr, "-E- %s line %d: Wrong number of dimensions for %s. Need 2 got %d.\n",
1506  __FILE__, __LINE__, sdsname, ndims);
1507  exit(1);
1508 
1509  }
1510  // printf("ANGLE_LIMIT: %s n_senz=%d n_solz=%d\n",infile,*n_senz,*n_solz);
1511 
1512  if ((angle_limit = (float *) calloc((*n_senz) * (*n_solz), sizeof (float))) == NULL) {
1513  printf("-E- : Error allocating memory to angle_limit\n");
1514  exit(FATAL_ERROR);
1515  }
1516 
1517  // for (i=0; i< *n_senz;i++) {
1518  // if ( (anglelimit[i] = (float *)calloc((*n_solz) ,sizeof(float))) == NULL) {
1519  // printf("-E- : Error allocating memory to tindx\n");
1520  // exit(FATAL_ERROR);
1521  // }
1522  // }
1523 
1524  if (nc_get_var(ncid, sds_id, angle_limit) != NC_NOERR) {
1525  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
1526  __FILE__, __LINE__, sdsname, infile);
1527  exit(1);
1528  }
1529 
1530  // for (j=0;j<*n_senz;j++)
1531  // for (i=0;i<*n_solz;i++) {
1532  // //anglelimit[j][i] = *(angle_limit+j*(*n_solz)+i);
1533  // printf("RJH: angle_limt: %f %f %f\n",solz[i],senz[j],*(angle_limit+j*(*n_solz)+i));
1534  //
1535  // }
1536 
1537  *insenz = (float *) senz;
1538  *insolz = (float *) solz;
1539  *insolz = (float *) solz;
1540  *anglelimit = (float *) angle_limit;
1541 
1542  } else {
1543  fprintf(stderr, "-E- %s line %d: Error opening infile = %s.\n",
1544  __FILE__, __LINE__, infile);
1545  return (1);
1546  }
1547 
1548 
1549  return (0);
1550 
1551 }
1552 
1553 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) {
1554  int i, j;
1555 
1556  i = *ii;
1557  j = *jj;
1558 
1559  if (i < 0 || i >= n_senz) i = 0;
1560  if (j < 0 || j >= n_solz) j = 0;
1561 
1562  if (insenz > senz[i]) {
1563  while (i < n_senz && insenz > senz[i])
1564  i++;
1565  if (i >= n_senz) i = n_senz - 1;
1566  } else {
1567  while (i >= 0 && insenz < senz[i])
1568  i--;
1569  if (i < 0) i = 0;
1570  }
1571 
1572  if (insolz > solz[j]) {
1573  while (j < n_solz && insolz > solz[j])
1574  j++;
1575  if (j >= n_solz) j = n_solz - 1;
1576  } else {
1577  while (j >= 0 && insolz < solz[j])
1578  j--;
1579  if (j < 0) j = 0;
1580  }
1581 
1582  *ii = i;
1583  *jj = j;
1584 
1585  return (anglelimit[i][j]);
1586 }
1587 
1588 /********************************************************************************
1589 !* *
1590 !* Name: FINDMATCH *
1591 !* Purpose: finds the closest match for ELEM in LIST *
1592 !* Parameters: LIST - array of values to match ELEM to. Elements is array *
1593 !* should increase in value. *
1594 !* Algorithm: linearly compare ELEM to each element. When ELEM is smaller *
1595 !* than the LIST(I), then it is assumed to be closest to LIST(I-1) *
1596 !* or LIST(I). The one that has the smallest absolute difference *
1597 !* to ELEM is returned as the closest match. *
1598 !* Globals used: none *
1599 !* Global output: none *
1600 !* Return Codes: the closest matching element index *
1601 !* Special Considerations: none *
1602 !* *
1603 !********************************************************************************/
1604 
1605 int32_t findMatch(float *list, int32_t nbands, float elem) {
1606  int i;
1607  float diff1, diff2;
1608 
1609  for (i = 0; i < nbands; i++)
1610  if (list[i] > elem) break;
1611 
1612  diff1 = fabs(list[i - 1] - elem);
1613  diff2 = fabs(list[i] - elem);
1614 
1615  if (diff1 < diff2)
1616  return (i - 1);
1617  else
1618  return (i);
1619 
1620 }
1621 
1622 /*
1623 !********************************************************************************
1624 !* *
1625 !* Name: channelRatio *
1626 !* Purpose: Calculate 3-channel ratios. *
1627 !* Parameters: none *
1628 !* Algorithm: The 0.94-um water vapor absorption channel is ratioed against *
1629 !* the linear combination of two window channels near 0.86 and *
1630 !* 1.03 um to obtain one channel ratio for the 0.94-um band. *
1631 !* Similar calculation is done for the 1.14-um water vapor band. *
1632 !* Globals used: NB1,NB2,NBP94,NB3,NB4,NB1P14 - number of points used in *
1633 !* channel ratios for both the .94- and 1.14-um regions*
1634 !* IST1,IED1,IST2,IED2,ISTP94,IEDP94 - 3-channel ratioing *
1635 !* parameters for the 0.94-um water vapor band *
1636 !* IST3,IED3,IST4,IED4,IST1P14,IED1P14 - 3-channel ratioing *
1637 !* parameters for the 1.14-um water vapor band. *
1638 !* WT1,WT2,WT3,WT4,JA - Relative weights for the four window *
1639 !* channels used in channel-ratioing calculations. JA *
1640 !* is an output parameter from a table searching *
1641 !* routine. *
1642 !* TRNCAL - Atmospheric transmittance spectra. *
1643 !* Global output:R094,R114 - 3-channel ratio values for the 0.94- and 1.14-um *
1644 !* water vapor bands. *
1645 !* Return Codes: none *
1646 !* Special Considerations: none *
1647 !* *
1648 !********************************************************************************
1649  */
1650 void channelRatio() {
1651 
1652  float const1[NH2OMAX], const2[NH2OMAX], const3[NH2OMAX], const4[NH2OMAX], const5[NH2OMAX], const6[NH2OMAX];
1653  int32_t nb1 = getinput7_.nb1, nb2 = getinput7_.nb2, nbp94 = getinput7_.nbp94, \
1654  nb3 = getinput7_.nb3, nb4 = getinput7_.nb4, nb1p14 = getinput7_.nb1p14;
1655  int32_t ist1 = init_speccal6_.ist1, ied1 = init_speccal6_.ied1, ist2 = init_speccal6_.ist2, \
1656  ied2 = init_speccal6_.ied2, istp94 = init_speccal6_.istp94, iedp94 = init_speccal6_.iedp94;
1657  int32_t ist3 = init_speccal7_.ist3, ied3 = init_speccal7_.ied3, ist4 = init_speccal7_.ist4, \
1658  ied4 = init_speccal7_.ied4, ist1p14 = init_speccal7_.ist1p14, ied1p14 = init_speccal7_.ied1p14;
1659  float wt1 = init_speccal8_.wt1, wt2 = init_speccal8_.wt2, wt3 = init_speccal8_.wt3, wt4 = init_speccal8_.wt4;
1660 
1661  int i, j;
1662 
1663  // Calculate average of spectra over window and water vapor absorption regions.
1664 
1665  for (j = 0; j < NH2OMAX; j++) {
1666  const1[j] = 0.0;
1667  const2[j] = 0.0;
1668  const3[j] = 0.0;
1669  const4[j] = 0.0;
1670  const5[j] = 0.0;
1671  const6[j] = 0.0;
1672 
1673  for (i = ist1 - 1; i < ied1; i++)
1674  const1[j] = const1[j] + tran_table1_.trntbl[j][i];
1675 
1676  const1[j] /= nb1;
1677 
1678  for (i = ist2 - 1; i < ied2; i++)
1679  const2[j] = const2[j] + tran_table1_.trntbl[j][i];
1680  const2[j] /= nb2;
1681 
1682  for (i = istp94 - 1; i < iedp94; i++)
1683  const3[j] = const3[j] + tran_table1_.trntbl[j][i];
1684  const3[j] /= nbp94;
1685 
1686  tran_table1_.r0p94[j] = const3[j] / (wt1 * const1[j] + wt2 * const2[j]);
1687 
1688  for (i = ist3 - 1; i < ied3; i++)
1689  const4[j] = const4[j] + tran_table1_.trntbl[j][i];
1690  const4[j] /= nb3;
1691 
1692  for (i = ist4 - 1; i < ied4; i++)
1693  const5[j] = const5[j] + tran_table1_.trntbl[j][i];
1694  const5[j] /= nb4;
1695 
1696  for (i = ist1p14 - 1; i < ied1p14; i++)
1697  const6[j] = const6[j] + tran_table1_.trntbl[j][i];
1698  const6[j] /= nb1p14;
1699 
1700  tran_table1_.r1p14[j] = const6[j] / (wt3 * const4[j] + wt4 * const5[j]);
1701  }
1702 
1703 
1704 }
1705 
1706 void kdistgasabs(float *kcdf, float *abscf, float*waveno, float *wavobs, int32_t nphi, int32_t n_layers, int32_t nbands) {
1707  //
1708  // R. Healy 3/2016 - converted from fortran kdist_gas_obs in atrem code which is based on Amir Ibrahim's
1709  // Matlab code
1710  // Construct k coefficients from gas transmission coefficients (from Amir Ibrahim's Matlab Code)
1711  //
1712  // abscf[nlayers][np_hi] (IN) :: the absorption coefficients
1713  // waveno[np_hi] (IN) :: wave number of the spectrum
1714  // wavobs[nwave] (IN) :: wavelengths of the instrument detector
1715  // np_hi (IN) :: Number of high resolution spectral points,
1716  // covering 3000 cm-1 to 18,000 cm-1 at
1717  // point spacing of 0.05 cm-1
1718  // nlayers (IN) :: number of atmospheric layers
1719  // nwave (IN) :: number of wavelengths of the instrument
1720  // kcdf[7][nwave][nlayers] (OUT):: K-Coefficients
1721 
1722  float **alayers, g7[7] = {0, 0.379, 0.6, 0.81, 0.9548, 0.9933, 1};
1723  float *diflam, *UV_lam, *IR_lam, *UV_dlam, *lam, *wn; //,*ir_dlam
1724  float *dwave, *dwn, *g_i, *k_i, *logk, *kint[7], **k7[7], *kk;
1725  int32_t **wavel_window, *binnum;
1726  float dlam, wmin, wmax, swav;
1727  float Q = {2.15199993E+25}; // Q=# of molecules above the surface at one atmosphere (molecules/cm**2)
1728  float AvN = {6.0225E+23};
1729  int kwavdn, kwavup, nbins;
1730  int32_t nsamp;
1731  int k, i, j, n, ndxwuv, ndxwir, ndxtot, iw;
1732 
1733  alayers = (float **) malloc(n_layers * sizeof (float *));
1734  for (i = 0; i < n_layers; i++) {
1735  if ((alayers[i] = (float *) malloc(nphi * sizeof (float))) == NULL) {
1736  printf("%s, %d - E - unable to allocate memory in kdistgasabs\n",
1737  __FILE__, __LINE__);
1738  exit(EXIT_FAILURE);
1739  }
1740  for (n = 0; n < nphi; n++)
1741  alayers[i][n] = *(abscf + i * n_layers + n) * Q * 28.966 / AvN / 1.0e-6;
1742  //alayers[i][n] = abscf[i][n]*Q*28.966/ AvN / 1.0e-6;
1743  }
1744 
1745  diflam = (float *) malloc(nbands * sizeof (float));
1746 
1747  dlam = 9999.;
1748  wmin = 9999;
1749  wmax = -9999;
1750  for (i = 0; i < nbands - 1; i++) {
1751  diflam[i] = wavobs[i + 1] - wavobs[i];
1752  if (dlam > diflam[i]) dlam = diflam[i];
1753  if (wmin > wavobs[i]) wmin = wavobs[i];
1754  if (wmax < wavobs[i]) wmax = wavobs[i];
1755  }
1756 
1757  diflam[nbands - 1] = diflam[nbands - 2];
1758 
1759  if (wmin > wavobs[nbands - 1]) wmin = wavobs[nbands - 1];
1760  if (wmax < wavobs[nbands - 1]) wmax = wavobs[nbands - 1];
1761 
1762  ndxwuv = (wmin - 0.0001 - 0.3) / dlam + 1;
1763  ndxwir = (3.1 - (wmax + 0.0001)) / diflam[nbands - 1] + 1;
1764  ndxtot = nbands + ndxwuv + ndxwir;
1765 
1766  UV_lam = (float *) malloc(ndxwuv * sizeof (float));
1767  IR_lam = (float *) malloc(ndxwir * sizeof (float));
1768  UV_dlam = (float *) malloc(ndxwuv * sizeof (float));
1769  lam = (float *) malloc(ndxtot * sizeof (float));
1770  dwave = (float *) malloc(ndxtot * sizeof (float));
1771  dwn = (float *) malloc(ndxtot * sizeof (float));
1772  wn = (float *) malloc(ndxtot * sizeof (float));
1773  wavel_window = (int **) malloc(2 * sizeof (int *));
1774  binnum = (int32_t *) malloc(n_layers * sizeof (int32_t));
1775 
1776  for (i = 0; i < 2; i++)
1777  wavel_window[i] = (int *) malloc(ndxtot * sizeof (int));
1778 
1779  k = ndxtot - 1;
1780  swav = 0.3;
1781 
1782  for (i = 0; i < ndxwuv; i++) {
1783  UV_lam[i] = swav;
1784  UV_dlam[i] = diflam[0];
1785  swav = swav + dlam;
1786  lam[k] = UV_lam[i];
1787  dwn[k] = 10000. * dlam / pow(lam[k], 2);
1788  wn[k] = 10000. / lam[k];
1789  k--;
1790  }
1791  swav = wmin - 1;
1792  for (i = 0; i < nbands; i++) {
1793  lam[k] = swav;
1794  dwn[k] = 10000. * diflam[i] / pow(lam[k], 2);
1795  wn[k] = 10000. / lam[k];
1796  swav = swav + diflam[i];
1797  k--;
1798  }
1799 
1800  swav = wmax + 0.0001;
1801 
1802  for (i = 0; i < ndxwir; i++) {
1803  IR_lam[i] = swav;
1804  //ir_dlam[i] = diflam[nbands];
1805  swav = swav + diflam[nbands - 1];
1806  lam[k] = IR_lam[i];
1807  dwn[k] = 10000. * diflam[nbands - 1] / pow(lam[k], 2);
1808  wn[k] = 10000. / lam[k];
1809  k--;
1810  }
1811 
1812  for (i = 0; i < 7; i++) {
1813  kint[i] = (float *) malloc(ndxtot * sizeof (float));
1814  k7[i] = (float **) malloc(ndxtot * sizeof (float *));
1815  for (k = 0; k < ndxtot; k++)
1816  k7[i][k] = (float *) malloc(n_layers * sizeof (float));
1817  }
1818 
1819  iw = nbands - 1; // index for wavelngths for output array kcdf
1820 
1821  for (i = ndxwir; i < ndxtot - ndxwuv; i++) {
1822  kwavdn = -999;
1823  kwavup = -999;
1824 
1825  for (k = 0; k < nphi; k++) {
1826  if (kwavdn < 0 && waveno[k]> (wn[i] - dwn[i] / 2.0)) kwavdn = k;
1827  if (kwavup < 0 && waveno[k]> (wn[i] + dwn[i] / 2.0)) kwavup = k;
1828  }
1829 
1830  wavel_window[0][i] = kwavdn;
1831  wavel_window[1][i] = kwavup;
1832 
1833  if (kwavup < 0 || kwavdn < 0)
1834  nsamp = 1;
1835  else
1836  nsamp = kwavup - kwavdn + 1;
1837 
1838  nbins = nsamp;
1839  g_i = (float *) malloc((nbins + 1) * sizeof (float));
1840  k_i = (float *) malloc((nbins + 1) * sizeof (float));
1841  kk = (float *) malloc((nsamp + 1) * sizeof (float));
1842  logk = (float *) malloc((nbins + 1) * sizeof (float));
1843 
1844  for (k = 0; k < n_layers; k++) {
1845  if (kwavdn < 0 || kwavup < 0)
1846  for (j = 0; j < nsamp + 1; j++)
1847  kk[j] = 0;
1848  else
1849  n = kwavdn - 1;
1850  for (j = 0; j < nsamp && n < nphi; j++) {
1851  kk[j] = alayers[k][n];
1852  n++;
1853  }
1854  binnum[k] = nbins;
1855  ecdf_(k_i, g_i, &binnum[k], kk, &nsamp);
1856 
1857  for (j = 0; j < nbins + 1; j++)
1858  logk[j] = log10(k_i[j]);
1859 
1860  for (j = 0; j < 7; j++) {
1861 
1862  kint[j][i] = linterp(g_i, logk, binnum[k], g7[j]);
1863  k7[j][i][k] = pow(10, kint[j][i]);
1864 
1865  if (i >= ndxwir && i < ndxtot - ndxwuv) {
1866  if (isnan(k7[j][i][k])) k7[j][i][k] = 0;
1867 
1868  //kcdf[j][iw][k] = k7[j][i][k];
1869  // *(kcdf+j*nbands*7+iw*n_layers+k) = k7[j][i][k];
1870  *(kcdf + j * nbands * n_layers + iw * n_layers + k) = k7[j][i][k];
1871  }
1872  }
1873 
1874  }
1875 
1876  free(logk);
1877  free(g_i);
1878  free(k_i);
1879  free(kk);
1880 
1881  if (i >= ndxwir && i < ndxtot - ndxwuv) iw--;
1882 
1883  }
1884 
1885  free(UV_lam);
1886  free(IR_lam);
1887  free(UV_dlam);
1888  //free ( ir_dlam );
1889  free(lam);
1890  free(dwave);
1891  free(dwn);
1892  free(wn);
1893  free(wavel_window[0]);
1894  free(wavel_window[1]);
1895  for (i = 0; i < 7; i++) {
1896  for (k = 0; k < ndxtot; k++)
1897  free(k7[i][k]);
1898  free(k7[i]);
1899  free(kint[i]);
1900  }
1901  for (i = 0; i < n_layers; i++)
1902  free(alayers[i]);
1903  free(alayers);
1904  free(binnum);
1905 }
1906 
1907 /*
1908  * Converted from Fortran by R. Healy (SAIC) 3/2016
1909 !********************************************************************************
1910 !* *
1911 !* Name: model_adjust *
1912 !* Purpose: resets the bottom boundary of the input model if the surface *
1913 !* elevation is greater than 0, and calculate the column water vapor *
1914 !* amount in the selected model. *
1915 !* Parameters: none *
1916 !* Algorithm: If the surface elevation > 0, the bottom layer temperature and *
1917 !* water vapor volume mixing ratio are obtained through linear *
1918 !* interpolation, while the bottom layer pressure is obtained *
1919 !* through exponential interpolation. *
1920 !* Globals used: H, T, P, VMR, NB, NL, HSURF - these values are adjusted if *
1921 !* HSURF > 0 *
1922 !* Global output: CLMVAP - Column water vapor amount in unit of cm. *
1923 !* Q - Number of molecules above the surface at one *
1924 !* atmosphere in units of molecules/cm**2 *
1925 !* Return Codes: none *
1926 !* Special Considerations: none *
1927 !* *
1928 !********************************************************************************
1929  */
1930 
1931 void model_adjust() {
1932 
1933  float *h = getinput3_.h, *t = getinput3_.t, *p = getinput3_.p, *vmr = getinput3_.vmr;
1934  int32_t nl = getinput3_.nl, nb = getinput3_.nb;
1935  float hsurf = getinput5_.hsurf;
1936  float xppp = getinput14_.xppp;
1937 
1938  static float hp[MODELMAX], tp[MODELMAX], pp[MODELMAX], vmrp[MODELMAX];
1939 
1940  float dhk, dhs, tsurf, vmrs, psurf, amtvrt, damtvrt, hplane;
1941  float q, tplane, pplane, vmrsp, dhss, dhkk, amtvrtp, damtvrtp;
1942  int32_t kk=-1, i, k_surf;
1943 
1944  q = 2.152e25;
1945  model_adj1_.q = q;
1946 
1947  //Determine index of H() such that H(index) < HSURF < H(index+1)
1948 
1949  for (i = 0; i < nb; i++)
1950  if (hsurf == h[i]) hsurf = h[i] + 0.0001;
1951 
1952  // Determine index of H() such that H(index) < HSURF < H(index+1)
1953  locate_pos_(h, &nb, &hsurf, &k_surf);
1954 
1955  // K_SURF is an index relative to the original model atmosphere (0 - 100 km)
1956  if (k_surf < 0) {
1957  fprintf(stderr, "***WARNING: Surface elevation smaller then lowest boundary of the model atmosphere.\n");
1958  k_surf = 0;
1959  } else {
1960  dhk = h[k_surf + 1] - h[k_surf];
1961  dhs = hsurf - h[k_surf];
1962  //linear interpolation for surface temperature (TSURF) and VMR ( VMRS)
1963  tsurf = t[k_surf]+(dhs / dhk)*(t[k_surf + 1] - t[k_surf]);
1964  vmrs = vmr[k_surf]+(dhs / dhk)*(vmr[k_surf + 1]);
1965  //exponential interpolation for surface pressure (PSURF)
1966  psurf = p[k_surf] * exp(-log(p[k_surf]));
1967  h[0] = hsurf;
1968  p[0] = psurf;
1969  t[0] = tsurf;
1970  vmr[0] = vmrs;
1971 
1972  nb -= k_surf;
1973  nl = nb - 1;
1974 
1975  for (i = 1; i < nb; i++) {
1976  h[i] = h[i + k_surf];
1977  p[i] = p[i + k_surf];
1978  t[i] = t[i + k_surf];
1979  vmr[i] = vmr[i + k_surf];
1980  }
1981  //Zero out pressures and VMRS of top atmospheric layers.
1982  for (i = nb; i < MODELMAX; i++) {
1983  h[i] = 1000;
1984  p[i] = 0.0;
1985  t[i] = 300.;
1986  vmr[i] = 0.0;
1987  }
1988  }
1989  amtvrt = 0.0;
1990 
1991  for (i = 0; i < nl; i++) {
1992  damtvrt = (q) *(p[i] - p[i + 1])*(vmr[i] + vmr[i + 1]) / 2.0;
1993  amtvrt += damtvrt;
1994  }
1995 
1996  model_adj1_.clmvap = amtvrt / 3.34e22;
1997 
1998  printf("ATREM: Model_adjust: Column vapor amount in model atmosphere from ground to space = %f cm\n", model_adj1_.clmvap);
1999  /*
2000  Setting the upward atmospheric path's T, P, and VMR profiles:
2001 
2002  1st duplicate the entire atmospheric profiles from the downward path
2003  to the upward path
2004 
2005  */
2006 
2007  for (i = 0; i < MODELMAX; i++) {
2008  hp[i] = h[i];
2009  tp[i] = t[i];
2010  pp[i] = p[i];
2011  vmrp[i] = vmr[i];
2012  }
2013  // Set the highest plane altitude to the upper bound of model atmosphere
2014  hplane = xppp;
2015 
2016  if (hplane >= 100.0) hplane = 100.0 - 0.0001;
2017  //Do special processing if the plane height (HPLANE) is greater than HP(1)
2018  if (hplane > hp[0]) {
2019  // Reset Plane altitude HPLANE (= XPPP) to a larger value if
2020  // HPLANE.EQ.HP(I) to avoid possible problems in table
2021  // searching using LOCATE
2022  for (i = 0; i < MODELMAX; i++)
2023  if (hplane == hp[i]) hplane = hp[i] - 0.0001;
2024  //Determine index of HP() such that HP(index) < HPLANE < H(index+1)
2025  locate_pos_(hp, &nb, &hplane, &kk);
2026 
2027  if (kk < 0)
2028  printf("WARNING: Plane altitude less then lowest boundary of the model atmosphere.\n");
2029  else {
2030  dhkk = hp[kk + 1] - hp[kk];
2031  dhss = hplane - hp[kk];
2032  //linear interpolation for plane temperature (TPLANE) and VMR ( VMRSP)
2033  tplane = tp[kk] + (dhss / dhkk)*(tp[kk + 1] - tp[kk]);
2034  vmrsp = vmrp[kk] + (dhss / dhkk)*(vmrp[kk + 1] - vmrp[kk]);
2035  //exponential interpolation for plane pressure (PPLANE)
2036  pplane = pp[kk] * exp(-log(pp[kk] / pp[kk + 1]) * dhss / dhkk);
2037  hp[kk + 1] = hplane;
2038  pp[kk + 1] = pplane;
2039  tp[kk + 1] = tplane;
2040  vmrp[kk + 1] = vmrsp;
2041  //Zero out pressures and VMRP of top atmospheric layers
2042  if (kk < MODELMAX - 2) {
2043  for (i = kk + 2; i < MODELMAX; i++) {
2044  hp[i] = 1000;
2045  pp[i] = 0.0;
2046  tp[i] = 300.;
2047  vmrp[i] = 0.0;
2048  }
2049  }
2050  }
2051  }
2052  amtvrtp = 0.0;
2053  for (i = 0; i < kk; i++) {
2054  damtvrtp = (q) *(pp[i] - pp[i + 1])*(vmrp[i] + vmrp[i + 1]) / 2.0;
2055  amtvrtp += damtvrtp;
2056  }
2057  model_adj3_.clmvapp = amtvrtp / 3.34e22;
2058  printf("ATREM: Model_adjust: Column vapor amount below plane = %f cm\n", model_adj3_.clmvapp);
2059 
2060  //Indices and parameters for the plane layer
2061  //model_adj3_.k_plane = kk;model_adj4_.k_surf = k_surf; // For final C version
2062  //printf("ATREM: MODEL_ADJUST - be sure to change next line to model_adj3_.k_plane = kk;model_adj4_.k_surf = k_surf once final conversion to C is complete\n ");
2063  model_adj3_.k_plane = kk + 1;
2064  model_adj4_.k_surf = k_surf + 1;
2065 
2066  model_adj3_.dvap_plane = q * (pp[kk] - pp[kk + 1]) * (vmrp[kk] + vmrp[kk + 1]) / 2.0 / 3.34e22;
2067  model_adj3_.dvap_layer = q * (p[kk] - p[kk + 1]) * (vmr[kk] + vmr[kk + 1]) / 2.0 / 3.34e22;
2068  model_adj3_.dp_plane = pp[kk] - pp[kk + 1];
2069  model_adj3_.dp_layer = p[kk] - p[kk + 1];
2070 
2071 }
2072 
2073 /*
2074  * ********************************************************************************
2075 !* *
2076 !* Name: LOCATE *
2077 !* Purpose: given an array XX of length N, and given a value X, returns a value*
2078 !* J such that X is between XX(J) and XX(J+1). XX must be monotonic, *
2079 !* Parameters: XX - monotonic array of values *
2080 !* N - number of elements in XX *
2081 !* X - value that will be matched to the XX array *
2082 !* J - index into the XX array where XX(J) <= X <= XX(J+1) *
2083 !* Algorithm: bisectional table searching, copied from Numerical Recipes. *
2084 !* Globals used: none *
2085 !* Global output: none *
2086 !* Return Codes: J=0 or J=N is returned to indicate that X is out of range *
2087 !* Special Considerations: none *
2088 !* *
2089 !********************************************************************************
2090  */
2091 void locate_pos_(float *xx, int32_t *n1, float *x1, int32_t *jj) {
2092  int32_t j, ju, jl, jm, n = *n1;
2093  float x = *x1;
2094 
2095  jl = -1;
2096  ju = n;
2097  while ((ju - jl) > 1) {
2098  jm = (ju + jl) / 2;
2099  if ((xx[n - 1] > xx[0]) && (x > xx[jm]))
2100  jl = jm;
2101  else
2102  ju = jm;
2103  }
2104  if (x == xx[0])
2105  j = -1;
2106  else if (x == xx[n - 1])
2107  j = n - 2;
2108  else
2109  j = jl;
2110 
2111  *jj = j;
2112 }
2113 
2114 /* Converted from Fortran by R. Healy 3/2016
2115  *
2116 !********************************************************************************
2117 !* *
2118 !* Name: GEOMETRY *
2119 !* Purpose: Calculates the solar and the observational geometric factors. *
2120 !* Parameters: none *
2121 !* Algorithm: The solar geometry was obtained based on the latitude, longitude,*
2122 !* GMT time using programs written by W. Mankin at National Center *
2123 !* for Atmospheric Research in Boulder, Colorado. The *
2124 !* geometric factors for CO2, O3, N2O, CO, CH4, and O2 are based *
2125 !* only on the solar and observational angles. Sixty artificial *
2126 !* geometric factors for H2O are set up to produce a transmittance *
2127 !* table for different atmospheric water vapor amounts. *
2128 !* Globals used: VRTO3 - Column O3 amount in units of atm-cm *
2129 !* IMN,IDY,IYR,IH,IM,IS - time and date of data measurements *
2130 !* XLATD,XLATM,XLATS,LATHEM - Latitude of measured area *
2131 !* XLONGD,XLONGM,XLONGS,LNGHEM - Longitude of measured area *
2132 !* CLMVAP - Column water vapor in unit of cm in the model atmosphere *
2133 !* Global output: *
2134 !* SOLZNI,SOLAZ,OBSZNI,OBSPHI,IDAY - Solar zenith angle, solar azimuth *
2135 !* angle, observational zenith angle, observational azimuth angle, *
2136 !* and the day number in the year *
2137 !* GCO2,GO3,GN2O,GCO,GCH4,GO2,SSH2O,GGEOM,TOTLO3 - Geometric factors for *
2138 !* different gases, and total O3 amount in the sun-surface ray path. *
2139 !* The geometric factor is defined as: if the vertical column amount *
2140 !* of the gas is equal 1, then GGAS is equal to the total amount of *
2141 !* the gas in the combined Sun-surface-sensor ray path. *
2142 !* Return Codes: none *
2143 !* Special Considerations: none *
2144 !* *
2145 !********************************************************************************
2146  */
2147 void geometry() {
2148  float vap_slant[NH2OMAX];
2149  int md[12] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}, leap_year;
2150  float ggeom;
2151  /*
2152  The VAPVRT array contains 60 column water vapor values used for generating
2153  a table of transmittance spectra with different amount of water vapor.
2154  The values in VAPVRT is designed so that there is approximately 2% change
2155  in the .94-um H2O channel ratio for column water vapor in the range .4-6 cm.
2156  */
2157  float vapvrt[NH2OMAX] = {.00, .02, .06, .11, .16, .21, .26, .31, .36, .40,
2158  .43, .46, .50, .54, .58, .62, .66, .70, .75, .80,
2159  .86, .92, .98, 1.06, 1.14, 1.22, 1.3, 1.4, 1.5, 1.6,
2160  1.7, 1.8, 1.9, 2.05, 2.2, 2.35, 2.55, 2.75, 2.95, 3.2,
2161  3.5, 3.8, 4.1, 4.4, 4.7, 5.0, 5.3, 5.6, 6.0, 6.4,
2162  7.0, 7.7, 8.5, 9.4, 10.4, 11.6, 13.0, 15.0, 25.0, 50.};
2163  float hplane = getinput14_.xppp;
2164  float obszni = geometry_l2gen_.senzn_l2, obsphi = geometry_l2gen_.senaz_l2, solzni = geometry_l2gen_.solzn_l2;
2165  float clmvap = model_adj1_.clmvap, clmvapp = model_adj3_.clmvapp;
2166  float *ssh2o = geometry2_.ssh2o;
2167  float dvap_plane = model_adj3_.dvap_plane, dp_plane = model_adj3_.dp_plane, dvap_layer = model_adj3_.dvap_layer, dp_layer = model_adj3_.dp_layer;
2168  float mu, mu0;
2169  int i, k_plane = model_adj3_.k_plane, splitpaths = geometry_l2gen_.splitpaths;
2170  float vap_sol, vap_sen;
2171  static int firstCall = 1;
2172 
2173  /*
2174 
2175  VAP_SLANT is a new array for containing two-way total vapor amounts
2176  Here the number "2" can be changed to other numbers, e.g., 2.5,
2177  without any major effects in derived water vapor values from
2178  imaging spectrometer data.
2179 
2180  */
2181 
2182  for (i = 0; i < NH2OMAX; i++)
2183  vap_slant[i] = vapvrt[i] * (float) 2.0;
2184 
2185  obszni /= RAD_DEG;
2186  obsphi /= RAD_DEG;
2187  solzni /= RAD_DEG;
2188 
2189  mu0 = (float) 1. / cos(solzni);
2190  mu = (float) 1. / cos(obszni);
2191  ggeom = mu0 + mu;
2192 
2193  geometry5_.mu0 = mu0;
2194  geometry5_.mu = mu;
2195  geometry2_.ggeom = ggeom;
2196  geometry1_.obszni = obszni;
2197  geometry1_.obsphi = obsphi;
2198  geometry1_.solzni = solzni;
2199 
2200  if (want_verbose || firstCall) printf("ATREM: GGEOM =%f OBSZNI = %f OBSPHI = %f solzni=%f degrees :: MU0=%f, MU = %f\n",
2201  ggeom, obszni, obsphi, solzni, mu0, mu);
2202  geometry2_.gco2 = ggeom;
2203 
2204  geometry2_.go3 = ggeom;
2205  if (hplane < 27.) geometry2_.go3 = ggeom - (float) 1. / cos(obszni);
2206 
2207  geometry2_.gn2o = ggeom;
2208  geometry2_.gco = ggeom;
2209  geometry2_.gch4 = ggeom;
2210  geometry2_.go2 = ggeom;
2211 
2212  geometry2_.totlo3 = geometry2_.go3 * getinput3_.vrto3;
2213 
2214  if (want_verbose || firstCall) printf("ATREM: TOTLO3 = %f %f\n", geometry2_.totlo3, getinput3_.vrto3);
2215 
2216  /*
2217  * Initialize newly created geometrical factors for each atmospheric
2218  layers (here G_VAP and G_OTHER are true geometrical factors that
2219  account for the actual Sun-surface-plane path lengths in the
2220  model atmosphere):
2221  *
2222  */
2223  //printf("ATREM: GEOMETRY - be sure to change line to k_plane NOT k_plane-1 once final conversion to C is complete\n ");
2224  //---For layers below the plane layer---
2225  //for (i=0;i<k_plane;i++) {
2226  for (i = 0; i < k_plane - 1; i++) {
2227  geometry3_.g_vap[i] = ggeom;
2228  geometry3_.g_other[i] = ggeom;
2229  }
2230 
2231  //printf("ATREM: GEOMETRY - be sure to change line to k_plane+1 NOT k_plane once final conversion to C is complete\n ");
2232  /*---For layers above the plane layer */
2233  // for (i=k_plane+1;i<MODELMAX;i++) {
2234  for (i = k_plane; i < MODELMAX; i++) {
2235  geometry3_.g_vap[i] = ggeom - (float) 1. / cos(obszni);
2236  geometry3_.g_other[i] = ggeom - (float) 1. / cos(obszni);
2237  }
2238  /* ---Special treatment for the plane layer to take account the
2239  "shorter" upward path length
2240  */
2241  //printf("ATREM: GEOMETRY - be sure to change line to k_plane NOT k_plane-1 once final conversion to C is complete\n ");
2242  /*
2243  geometry3_.g_vap[k_plane] = ggeom - 1./cos(obszni) \
2244  + dvap_plane/dvap_layer/cos(obszni);
2245  geometry3_.g_other[k_plane] = ggeom - 1./cos(obszni) \
2246  + dp_plane/dp_layer/cos(obszni);
2247  */
2248  geometry3_.g_vap[k_plane - 1] = ggeom - 1. / cos(obszni)
2249  + dvap_plane / dvap_layer / cos(obszni);
2250  geometry3_.g_other[k_plane - 1] = ggeom - 1. / cos(obszni)
2251  + dp_plane / dp_layer / cos(obszni);
2252  /*
2253  Calculate the water vapor SCALING factor relative to the total amount
2254  of water vapor in the model atmosphere in the L-shaped
2255  Sun-surface-plane ray path.
2256  */
2257  geometry4_.vap_slant_mdl = clmvap / cos(solzni) + clmvapp / cos(obszni);
2258  vap_sol = clmvapp*mu0;
2259  vap_sen = clmvapp*mu;
2260  if (want_verbose || firstCall) printf("ATREM: VAP_SLANT_MDL = %f cm\n", geometry4_.vap_slant_mdl);
2261 
2262  /*
2263  The "equivalent" geometrical factor corresponding to the total
2264  slant vapor amount of VAP_SLANT_MDL':
2265  */
2266 
2267  geometry3_.g_vap_equiv = geometry4_.vap_slant_mdl / clmvap;
2268 
2269  if (want_verbose || firstCall)
2270  printf("ATREM: G_VAP_EQUIV = %f clmvap=%f VAP_SOL=%f,VAP_SEN=%f \n", geometry3_.g_vap_equiv, clmvap,
2271  vap_sol, vap_sen);
2272 
2273  for (i = 0; i < NH2OMAX; i++) {
2274  ssh2o[i] = vap_slant[i] / geometry4_.vap_slant_mdl;
2275  // if (geometry_l2gen_.water_vapor > 0) {
2276  if (splitpaths) {
2277  geometry5_.ssh2o_s[0][i] = vapvrt[i] / vap_sol;
2278  geometry5_.ssh2o_s[1][i] = vapvrt[i] / vap_sen;
2279  // printf("SSH2O_S(1): %d %g %g | %g | %g \n",i,geometry5_.ssh2o_s[0][i],vapvrt[i],vap_sol,solzni*RAD_DEG);
2280  // printf("SSH2O_S(2): %d %g %g | %g | %g \n",i,geometry5_.ssh2o_s[1][i],vapvrt[i],vap_sen,obszni*RAD_DEG);
2281  }
2282  // } else if (splitpaths != 0) {
2283  // printf("ATREM: Split paths is not working because WaterVapor is 0\n");
2284  // }
2285  }
2286  /*
2287  Calculate the number of days that have passed in this year. Take leap year
2288  into account.
2289  */
2290  geometry1_.day = md[getinput8_.imn] + getinput8_.idy;
2291  leap_year = getinput8_.iyr - (4 * (getinput8_.iyr / 4));
2292  if ((!leap_year) && (geometry1_.day > 59) && (getinput8_.imn != 2)) geometry1_.day++;
2293  firstCall = 0;
2294 }
2295 
2296 void get_abscf_data(int levels, int bands, int sds_id, char filename[FILENAME_MAX], float* abscf, char *varname) {
2297  size_t start[3], count[3];
2298  int retval;
2299  int xid;
2300 
2301  start[0] = 0;
2302  start[1] = 0;
2303  start[2] = 0;
2304  count[0] = levels;
2305  count[1] = bands;
2306  count[2] = 0;
2307  retval = nc_inq_varid(sds_id, varname, &xid);
2308  if (retval != NC_NOERR) {
2309  fprintf(stderr,
2310  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
2311  __FILE__, __LINE__, filename, varname);
2312  exit(FATAL_ERROR);
2313  }
2314  retval = nc_get_vara_float(sds_id, xid, start, count, abscf);
2315  if (retval != NC_NOERR) {
2316  fprintf(stderr,
2317  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
2318  __FILE__, __LINE__, filename, varname);
2319  exit(FATAL_ERROR);
2320  }
2321 
2322  return;
2323 }
2324 
2325 void apply_gas_trans(int32_t k_surf, int levels, int bands, float* abscf, float* dp, float* vmrm, float *tran_hi) {
2326  static float *sumcf;
2327  int32_t i, j;
2328  float Q = 2.152e25;
2329 
2330  if ((sumcf = (float *) calloc(bands, sizeof (float))) == NULL) {
2331  printf("-E- : Error allocating memory to sumcf\n");
2332  exit(FATAL_ERROR);
2333  }
2334 
2335  //printf("ATREM: INIT_SPECCAL - be sure to change below line to i=k_surf and i - k_surf (NOT i - k_surf + 1) once final conversion to C is complete\n ");
2336 
2337  for (i = 0; i < bands; i++)
2338  sumcf[i] = 0.;
2339 
2340  for (i = k_surf - 1; i < levels; i++)
2341  for (j = 0; j < bands; j++)
2342  sumcf[j] -= abscf[bands * i + j] * dp[i - k_surf + 1] * vmrm[i - k_surf + 1];
2343 
2344  for (i = 0; i < bands; i++)
2345  tran_hi[i] *= exp(sumcf[i] * Q * (float) 28.966 / (float) 6.0225E+23 / (float) 1.0E-06);
2346 
2347  free(sumcf);
2348 }
2349 
2350 /* Ported from Fortran code by R. Healy 3/17/2016
2351  * Comments from Fortran Code - unless noted otherwise
2352 !********************************************************************************
2353 !* *
2354 !* Name: INIT_SPECCAL *
2355 !* Purpose: initialize global data for spectrum calculations. *
2356 !* Parameters: none *
2357 !* Algorithm: initialize data. *
2358 !* Globals used: AH2O, APH2O, BH2O, BPH2O, SODLT, SOGAM, O3CF - Band model *
2359 !* parameters for spectral calculations. *
2360 !* WNDOW1, WNDOW2, WP94C, WNDOW3, WNDOW4, W1P14C - Center *
2361 !* positions of window and water vapor absorption *
2362 !* channels used in 3-channel ratio calculations. *
2363 !* NB1, NB2, NBP94, NB3, NB4, NB1P14 - Number of narrow channels *
2364 !* to form broader window and absorption channels. *
2365 !* Global output: *
2366 !* IH2OLQ,RLQAMT,NGASTT,NH2O,VSTART,VEND - Flag for including liquid *
2367 !* water, liquid water amount (cm), total number of gases *
2368 !* (typically 8), number of water vapor values, starting *
2369 !* and ending wavelengths in internal calculations. *
2370 !* NO3PT,NCV,NCVHAF,NCVTOT,VMIN,ISTART,IEND - Number of O3 abs. coef. *
2371 !* points, parameters for gaussian function and spectral *
2372 !* calculations. *
2373 !* ISTCAL,IEDCAL,DP,PM,TM,VMRM - Parameters for spectral calculations *
2374 !* IST1,IED1,IST2,IED2,ISTP94,IEDP94 - 3-channel ratioing parameters for *
2375 !* the 0.94-um water vapor band. *
2376 !* IST3,IED3,IST4,IED4,IST1P14,IED1P14 - 3-channel ratioing parameters for *
2377 !* the 1.14-um water vapor band. *
2378 !* WT1,WT2,WT3,WT4,JA - Relative weights for the four window channels *
2379 !* used in channel-ratioing calculations. JA is a *
2380 !* output parameter from a table searching routine. *
2381 !* NCV2,NCVHF2,NCVTT2,ISTRT2,IEND2,FINST2 - Parameters for smoothing *
2382 !* output reflectance spectra. *
2383 !* NATOT,NBTOT,NCTOT,NDTOT - Number of channels for the four AVIRIS' *
2384 !* grating spectrometers (A, B, C, and D). *
2385 !* Return Codes: None. *
2386 !* Special Considerations: Some parameters may need to be fine-tuned. *
2387 !* *
2388 !********************************************************************************
2389 !* Notes about water vapor VMRS and related quantities: *
2390 !* *
2391 !* VAPVRT(60) - a table containing 60 column vapor values (in unit of cm) *
2392 !* *
2393 !* VAP_SLANT(I) = VAPVRT(I) * 2.0, VAP_SLANT is a new table for containing *
2394 !* two-way total vapor amounts. Here the number "2" can be *
2395 !* changed to other numbers, e.g., 2.5, without major *
2396 !* effects on retrieved water vapor values. *
2397 !* *
2398 !* G_VAP(I = 1,..., NL) = true vapor geometric factor for each layer in *
2399 !* the model atmosphere (after adjusting for the elevated *
2400 !* surface. *
2401 !* *
2402 !* VMRM(I) = VMRM(I)*G_VAP(I). The VMRS are multiplied by the geometrical *
2403 !* factor. We can calculate the vapor transmittance on the *
2404 !* Sun-surface-sensor path by assuming a vertical path in *
2405 !* the model atmosphere with geometric-factor-adjusted VMRS. *
2406 !* *
2407 !* CLMVAP = vertical column amount from ground to space in model atmosphere*
2408 !* CLMVAPP = vertical column amount from ground to aircraft or satellite *
2409 !* sensor in model atmosphere *
2410 !* Q = 2.152E25 = # of molecules above the surface at one atmosphere *
2411 !* (in unit of molecules/cm**2) *
2412 !* *
2413 !* VAP_SLANT_MDL= CLMVAP/COS(SOLZNI) + CLMVAPP/COS(OBSZNI) = total amount *
2414 !* of water vapor in the model atmosphere in the L-shaped *
2415 !* Sun-surface-plane ray path. *
2416 !* *
2417 !* G_VAP_EQUIV = VAP_SLANT_MDL / CLMVAP = the "equivalent" geometrical *
2418 !* factor corresponding to the total slant vapor amount *
2419 !* VAP_SLANT_MDL and the column vapor amount CLMVAP. *
2420 !* *
2421 !* SSH2O(I) (I = 1, ..., 60) - a pure scaling factor relative to the total *
2422 !* slant vapor amount of VAP_SLANT_MDL, and *
2423 !* SSH2O(I) = VAP_SLANT(I) / VAP_SLANT_MDL *
2424 !* *
2425 !* SH2O = one value of SSH2O(I). SH2O is used during generation of the *
2426 !* look-up table. *
2427 !* *
2428 !* VAPTT = VAP_SLANT_MDL*SH2O, is the absolute total vapor amount on the *
2429 !* L-shaped path corresponding to a spectrum stored in the *
2430 !* look-up table. *
2431 !* *
2432 !* CLMWVP = 0.5*(VAPTTA+VAPTTB)/G_VAP_EQUIV, is the retrieved column water *
2433 !* vapor amount from imaging spectrometer data. *
2434 !********************************************************************************
2435  */
2436 void init_spectral_calculations() {
2437  int32_t co2 = getinput1_.co2, o2 = getinput1_.o2, n2o = getinput1_.n2o, co = getinput1_.co, ch4 = getinput1_.ch4, no2 = getinput1_.no2, o3 = getinput1_.o3;
2438  float *fwhm = getinput4_.fwhm, *wavobs = getinput4_.wavobs, *dp = init_speccal5_.dp, *tm = init_speccal5_.tm, *pm = init_speccal5_.pm, *vmrm = init_speccal5_.vmrm;
2439  float *finst2 = init_speccal10_.finst2;
2440  int32_t nbands = getinput5_.nbands;
2441  ;
2442  float totlo3 = geometry2_.totlo3; //=no2cf_init1_.rno2cf;
2443  float dlt2 = getinput5_.dlt2, wavcv2, dwvavr, sumins;
2444  float sno2 = getinput3_.sno2, totno2;
2445  float go3 = geometry2_.go3, Q = model_adj1_.q;
2446  static float *abscf_co2, *abscf_o2, *abscf_n2o, *abscf_ch4, *abscf_co, const1, *o3cf, *rno2cf; //=o3cf_init1_.o3cf;
2447  static int firstCall = 1;
2448  int32_t i, j, k_surf = model_adj4_.k_surf, nl = getinput3_.nl;
2449  int32_t iwndw1, iwndw2, iwndw4, iwndw5, iwp94c, iw1p14c, nchnla, nchnlb, nchnlc, nchnld, nb1haf, nb2haf, nb3haf, nb4haf;
2450  char *abscf_file = "abscf_gas.nc";
2451  char *abscf_no2o3file = "abscf_no2o3.nc";
2452  static size_t ncf;
2453 
2454  char *filedir;
2455  char filename[FILENAME_MAX];
2456  size_t start[3], count[3];
2457  int xid, yid;
2458  int retval, sds_id;
2459  static size_t bands, levels;
2460 
2461  float vrtno2 = 5.0E+15, gno2, sclgas;
2462  float *sumcf;
2463 
2464  if (firstCall) {
2465  const1 = 4.0 * log(2); // gaussian curve area constant
2466  init_speccal3_.nh2o = NH2OMAX; // number of water vapor values in table
2467  //Wavelength of medium resolution spectrum FWHM=.2 nm, .56-3.1 um
2468 
2469  for (i = 0; i < NP_MED; i++) {
2470  init_speccal12_.wavln_med[i] = (float) VSTART + i * (float) DWAVLN;
2471  /*
2472  Note: The grids of WAVNO_HI do not match the grids of 10000./WAVLN_MED.
2473  INDEX_MED is a specially designed index for finding close matches
2474  between the two kinds of grids.
2475  */
2476  init_speccal13_.index_med[i] = ((float) 10000. / init_speccal12_.wavln_med[i] - (float) 3000.) / (float) DWAVNO + 1;
2477  /*
2478  Note: WAVLN_MED_INDEX(I) is very close to WAVLN_MED(I),
2479  and WAVLN_MED_INDEX(I) >= WAVLN_MED(I)
2480 
2481 
2482  */
2483  init_speccal13_.wavln_med_index[i] = (float) 10000. / ((init_speccal13_.index_med[i] - (float) 1.)*(float) DWAVNO + (float) 3000.);
2484  }
2485  //Wavelength of medium resolution spectrum FWHM=.2 nm, .3-3.1 um
2486  for (i = 0; i < NP_STD; i++)
2487  init_speccal12_.wavln_std[i] = (float) 0.3 + i * (float) DWAVLN;
2488  /*
2489  Initialize arrays for smoothing medium resolution spectrum (DLT_MED = 0.2 nm,
2490  and point spacing DWAVLN = 0.0001 micron) to coarser spectral
2491  resolution data from imaging spectrometers.
2492 
2493  NCVHF(I) = ( FACDLT * FWHM(I) / DWAVLN + 1.)
2494 
2495  */
2496  for (i = 0; i < nbands; i++)
2497  init_speccal15_.ncvhf[i] = (float) FACDLT * fwhm[i] / (float) DWAVLN + (float) 1.0;
2498 
2499  // parameters and arrays to smooth output surface reflectance spectrum
2500  wavcv2 = FACDLT*dlt2;
2501  /*
2502  * Find the largest value in the FWHM array, and use this value in calculation
2503  of indices for smoothing output reflectance spectra. This smoothing
2504  algorithm should work well with grating spectrometers having nearly
2505  constant spectral resolutions, but not so well for prism spectrometers
2506  having variable spectral resolution.
2507  */
2508  dwvavr = fwhm[0];
2509 
2510  for (i = 1; i < nbands; i++)
2511  if (dwvavr < fwhm[i]) dwvavr = fwhm[i];
2512 
2513  init_speccal10_.ncv2 = wavcv2 / dwvavr;
2514  init_speccal10_.ncvhf2 = init_speccal10_.ncv2 + 1;
2515  init_speccal10_.ncvtt2 = init_speccal10_.ncv2 + 1;
2516 
2517  if (dlt2 != 0) {
2518  sumins = 0;
2519  for (i = init_speccal10_.ncvhf2 - 1; i < init_speccal10_.ncvtt2; i++) {
2520  finst2[i] = exp(-const1 * pow((i - init_speccal10_.ncvhf2 + 1) * dwvavr / dlt2, 2));
2521  sumins += finst2[i];
2522  }
2523  for (i = 0; i < init_speccal10_.ncvhf2 - 1; i++) {
2524  finst2[i] = finst2[init_speccal10_.ncvtt2 - i - 1];
2525  sumins += finst2[i];
2526  }
2527  sumins *= dwvavr;
2528 
2529  for (i = 0; i < init_speccal10_.ncvtt2; i++) {
2530  finst2[i] *= (dwvavr / sumins);
2531  }
2532  }
2533 
2534  init_speccal10_.istrt2 = init_speccal10_.ncvhf2;
2535  init_speccal10_.iend2 = nbands - init_speccal10_.ncvhf2;
2536 
2537  /* number of channels of the four AVIRIS spectrometers. These are used
2538  in removing null AVIRIS radiance values in the overlap portions of two
2539  adjacent spectrometers.
2540 
2541  Note by R.Healy - The Aviris correction isn't being used in l2gen yet. (3/2016)
2542  There are no sensor specific modifications being applied. Atrem has only been
2543  tested on HICO data.
2544  */
2545  nchnla = 32;
2546  nchnlb = 64;
2547  nchnlc = 64;
2548  nchnld = 64;
2549 
2550  init_speccal11_.natot = nchnla;
2551  init_speccal11_.nbtot = nchnla + nchnlb;
2552  init_speccal11_.nctot = nchnla + nchnlb + nchnlc;
2553  init_speccal11_.ndtot = nchnla + nchnlb + nchnlc + nchnld;
2554 
2555  /* Resetting window wavelength positions and calculating weights for
2556  window and absorption channels used in 3-channel ratioing.
2557  Note that the C version of this function returns 0 based indices
2558  */
2559  iwndw1 = findMatch(wavobs, nbands, getinput6_.wndow1);
2560  iwndw2 = findMatch(wavobs, nbands, getinput6_.wndow2);
2561 
2562  getinput6_.wndow1 = wavobs[iwndw1];
2563  getinput6_.wndow2 = wavobs[iwndw2];
2564 
2565  if ((getinput7_.nb1 % 2) == 0) getinput7_.nb1 += 1;
2566  if ((getinput7_.nb2 % 2) == 0) getinput7_.nb2 += 1;
2567 
2568  nb1haf = (getinput7_.nb1 - 1) / 2;
2569  nb2haf = (getinput7_.nb2 - 1) / 2;
2570 
2571  // Adding one to all the indexes because of legacy fortran - subtracted in get_atrem_cor
2572  init_speccal6_.ist1 = iwndw1 - nb1haf + 1;
2573  init_speccal6_.ied1 = iwndw1 + nb1haf + 1;
2574  init_speccal6_.ist2 = iwndw2 - nb2haf + 1;
2575  init_speccal6_.ied2 = iwndw2 + nb2haf + 1;
2576 
2577  iwp94c = findMatch(wavobs, nbands, getinput6_.wp94c);
2578  getinput6_.wp94c = wavobs[iwp94c];
2579 
2580  if ((getinput7_.nbp94 % 2) == 0) getinput7_.nbp94 += 1;
2581 
2582  nb3haf = (getinput7_.nbp94 - 1) / 2;
2583  init_speccal6_.istp94 = iwp94c - nb3haf + 1;
2584  init_speccal6_.iedp94 = iwp94c + nb3haf + 1;
2585 
2586  init_speccal8_.wt1 = (getinput6_.wndow2 - getinput6_.wp94c) / (getinput6_.wndow2 - getinput6_.wndow1);
2587  init_speccal8_.wt2 = (getinput6_.wp94c - getinput6_.wndow1) / (getinput6_.wndow2 - getinput6_.wndow1);
2588 
2589  iwndw4 = findMatch(wavobs, nbands, getinput6_.wndow3);
2590  iwndw5 = findMatch(wavobs, nbands, getinput6_.wndow4);
2591 
2592  getinput6_.wndow3 = wavobs[iwndw4];
2593  getinput6_.wndow4 = wavobs[iwndw5];
2594 
2595  if ((getinput7_.nb3 % 2) == 0) getinput7_.nb3 += 1;
2596  if ((getinput7_.nb4 % 2) == 0) getinput7_.nb4 += 1;
2597 
2598  nb3haf = (getinput7_.nb3 - 1) / 2;
2599  nb4haf = (getinput7_.nb4 - 1) / 2;
2600 
2601  init_speccal7_.ist3 = iwndw4 - nb3haf + 1;
2602  init_speccal7_.ied3 = iwndw4 + nb3haf + 1;
2603  init_speccal7_.ist4 = iwndw5 - nb4haf + 1;
2604  init_speccal7_.ied4 = iwndw5 + nb4haf + 1;
2605 
2606  iw1p14c = findMatch(wavobs, nbands, getinput6_.w1p14c);
2607  getinput6_.w1p14c = wavobs[iw1p14c];
2608 
2609  if ((getinput7_.nb1p14 % 2) == 0) getinput7_.nb1p14 += 1;
2610 
2611  nb3haf = (getinput7_.nb1p14 - 1) / 2;
2612  init_speccal7_.ist1p14 = iw1p14c - nb3haf + 1;
2613  init_speccal7_.ied1p14 = iw1p14c + nb3haf + 1;
2614 
2615  init_speccal8_.wt3 = (getinput6_.wndow4 - getinput6_.w1p14c) / (getinput6_.wndow4 - getinput6_.wndow3);
2616  init_speccal8_.wt4 = (getinput6_.w1p14c - getinput6_.wndow3) / (getinput6_.wndow4 - getinput6_.wndow3);
2617 
2618  // Open the netcdf4 input file
2619  if ((filedir = getenv("OCDATAROOT")) == NULL) {
2620  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
2621  exit(FATAL_ERROR);
2622  }
2623 
2624  strcpy(filename, filedir);
2625  strcat(filename, "/common/");
2626  strcat(filename, abscf_file);
2627  retval = nc_open(filename, NC_NOWRITE, &sds_id);
2628  if (retval != NC_NOERR) {
2629  fprintf(stderr, "-E- %s line %d: nc_open(%s) failed.\n",
2630  __FILE__, __LINE__, filename);
2631  exit(FATAL_ERROR);
2632  }
2633  // Now read in the absorption coefficients
2634  retval = nc_inq_dimid(sds_id, "level", &xid)
2635  || nc_inq_dimlen(sds_id, xid, &levels)
2636  || nc_inq_dimid(sds_id, "band", &yid)
2637  || nc_inq_dimlen(sds_id, yid, &bands);
2638 
2639  if (retval) {
2640  fprintf(stderr, "-E- %s line %d: nc_inq_dimid(%s) failed.%d %d %d \n",
2641  __FILE__, __LINE__, filename, (int) levels, (int) bands, retval);
2642  exit(FATAL_ERROR);
2643  }
2644 
2645  abscf_o2 = (float *) calloc(levels*bands, sizeof (float));
2646  if (!(abscf_o2)) {
2647  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for abscf_n2o.\n",
2648  __FILE__, __LINE__);
2649  exit(FATAL_ERROR);
2650  }
2651 
2652  abscf_co2 = (float *) calloc(levels*bands, sizeof (float));
2653  if (!(abscf_co2)) {
2654  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for abscf_co2.\n",
2655  __FILE__, __LINE__);
2656  exit(FATAL_ERROR);
2657  }
2658  abscf_co = (float *) calloc(levels*bands, sizeof (float));
2659  if (!(abscf_co)) {
2660  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for abscf_co.\n",
2661  __FILE__, __LINE__);
2662  exit(FATAL_ERROR);
2663  }
2664  abscf_ch4 = (float *) calloc(levels*bands, sizeof (float));
2665  if (!(abscf_ch4)) {
2666  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for abscf_ch4.\n",
2667  __FILE__, __LINE__);
2668  exit(FATAL_ERROR);
2669  }
2670  abscf_n2o = (float *) calloc(levels*bands, sizeof (float));
2671  if (!(abscf_n2o)) {
2672  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for abscf_n2o.\n",
2673  __FILE__, __LINE__);
2674  exit(FATAL_ERROR);
2675  }
2676 
2677  get_abscf_data(levels, bands, sds_id, filename, abscf_o2, "abscf_o2");
2678  get_abscf_data(levels, bands, sds_id, filename, abscf_co2, "abscf_co2");
2679  get_abscf_data(levels, bands, sds_id, filename, abscf_co, "abscf_co");
2680  get_abscf_data(levels, bands, sds_id, filename, abscf_ch4, "abscf_ch4");
2681  get_abscf_data(levels, bands, sds_id, filename, abscf_n2o, "abscf_n2o");
2682 
2683  nc_close(sds_id);
2684 
2685  strcpy(filename, filedir);
2686  strcat(filename, "/common/");
2687  strcat(filename, abscf_no2o3file);
2688  strcat(filename, "\0");
2689  retval = nc_open(filename, NC_NOWRITE, &sds_id);
2690  if (retval != NC_NOERR) {
2691  fprintf(stderr, "-E- %s line %d: nc_open(%s) failed.\n",
2692  __FILE__, __LINE__, filename);
2693  exit(FATAL_ERROR);
2694  }
2695  // Now read in the absorption coefficients
2696  retval = nc_inq_dimid(sds_id, "ncf", &xid)
2697  || nc_inq_dimlen(sds_id, xid, &ncf);
2698 
2699  if (retval) {
2700  fprintf(stderr, "-E- %s line %d: nc_inq_dimid(%s) failed. %d %d \n",
2701  __FILE__, __LINE__, filename, (int) ncf, retval);
2702  exit(FATAL_ERROR);
2703  }
2704 
2705  o3cf = (float *) calloc(ncf, sizeof (float));
2706  if (!(o3cf)) {
2707  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for abscf_o3.\n",
2708  __FILE__, __LINE__);
2709  exit(FATAL_ERROR);
2710  }
2711  rno2cf = (float *) calloc(ncf, sizeof (float));
2712  if (!(o3cf)) {
2713  fprintf(stderr, "-E- %s line %d: Failed to allocate memory for abscf_no2.\n",
2714  __FILE__, __LINE__);
2715  exit(FATAL_ERROR);
2716  }
2717 
2718  start[0] = 0;
2719  start[1] = 0;
2720  start[2] = 0;
2721  count[0] = ncf;
2722  count[1] = 0;
2723  count[2] = 0;
2724 
2725  retval = nc_inq_varid(sds_id, "o3cf", &xid);
2726 
2727  if (retval != NC_NOERR) {
2728  fprintf(stderr,
2729  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
2730  __FILE__, __LINE__, filename, "o3cf");
2731  exit(FATAL_ERROR);
2732  }
2733  retval = nc_get_vara_float(sds_id, xid, start, count, o3cf);
2734  if (retval != NC_NOERR) {
2735  fprintf(stderr,
2736  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
2737  __FILE__, __LINE__, filename, "o3cf");
2738  exit(FATAL_ERROR);
2739  }
2740 
2741  retval = nc_inq_varid(sds_id, "no2cf", &xid);
2742 
2743  if (retval != NC_NOERR) {
2744  fprintf(stderr,
2745  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
2746  __FILE__, __LINE__, filename, "no2cf");
2747  exit(FATAL_ERROR);
2748  }
2749  retval = nc_get_vara_float(sds_id, xid, start, count, rno2cf);
2750  if (retval != NC_NOERR) {
2751  fprintf(stderr,
2752  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
2753  __FILE__, __LINE__, filename, "no2cf");
2754  exit(FATAL_ERROR);
2755  }
2756 
2757 
2758  nc_close(sds_id);
2759 
2760  firstCall = 0;
2761  }
2762 
2763  // Initialize the TRAN_HI array for high resolution spectrum:
2764  for (i = 0; i < bands; i++)
2765  init_speccal1_.tran_hi_others[i] = 1.0;
2766 
2767  /*
2768  Calculate medium resolution O3 transmittances (0.2 nm, 2-way) with
2769  a point spacing of 0.1 nm between 0.3 and 0.8 micron.
2770  */
2771  if (o3) {
2772  for (i = 0; i < NO3PT; i++) {
2773  init_speccal16_.tran_o3_std[i] = exp(-totlo3 * o3cf[i]);
2774  //printf("TRAN_O3_STD %d %g %f %f\n",i+1,init_speccal16_.tran_o3_std[i],totlo3,o3cf[i]);
2775  }
2776  } else {
2777  /* If NO2 is not intended to be included in total atmospheric gaseous
2778  transmittance calculations, assigning TRAN_NO2_STD = 1.0:
2779  */
2780  for (i = 0; i < NO3PT; i++)
2781  init_speccal16_.tran_o3_std[i] = 1.0;
2782 
2783  }
2784  /*
2785  Calculate medium resolution NO2 transmittances (0.2 nm, 2-way) with
2786  a point spacing of 0.1 nm between 0.3 and 0.8 micron.
2787  */
2788 
2789  vrtno2 = 5.0E+15;
2790  vrtno2 *= sno2;
2791 
2792  gno2 = go3;
2793  totno2 = gno2 * vrtno2;
2794 
2795  if (no2) {
2796  for (i = 0; i < NO2PT; i++) {
2797  init_speccal17_.tran_no2_std[i] = exp(-totno2 * rno2cf[i]);
2798  //printf("TRAN_NO2_STD %d %g %g %g\n",i+1,init_speccal17_.tran_no2_std[i],totno2,rno2cf[i]);
2799  }
2800  } else {
2801  /* If NO2 is not intended to be included in total atmospheric gaseous
2802  transmittance calculations, assigning TRAN_NO2_STD = 1.0:
2803  */
2804  for (i = 0; i < NO2PT; i++)
2805  init_speccal17_.tran_no2_std[i] = 1.0;
2806  }
2807 
2808  for (i = 0; i < nl; i++) {
2809  dp[i] = getinput3_.p[i] - getinput3_.p[i + 1];
2810  pm[i] = (getinput3_.p[i] + getinput3_.p[i + 1]) / (float) 2.0;
2811  tm[i] = (getinput3_.t[i] + getinput3_.t[i + 1]) / (float) 2.0;
2812  }
2813 
2814  /*
2815 
2816 !C
2817 !C Calculate high resolution transmittances (0.05 cm-1) of CO2, N2O, CO,
2818 !C CH4, and O2 in the 0.56 - 3.1 micron range, and save values for
2819 !C calculating total atmospheric transmittances later.
2820 !C Because water vapor amounts are allowed to vary,
2821 !C the high resolution water vapor transmittances are calculated
2822 !C in subroutines TRAN_TABLE and TRANCAL. TRAN_TABLE provides variable
2823 !C water vapor amounts, and calls TRANCAL for the calculation of
2824 !C corresponding vapor transmittance spectrum.
2825 !C
2826 
2827  */
2828 
2829  // For CO2 transmittance calculation -
2830  if (co2) {
2831  /*
2832  -- On 2/7/2013 B.-C. Gao made the modification - Increased SCLCO2 i
2833  -- from 1.0 to 1.1 to reflect the fact that the CO2 VMR reached the
2834  -- 2012 level of 390 ppmv.
2835  --
2836  */
2837  sclgas = 1.1;
2838  /*
2839  Scale the VMRM by the two-way path geometrical factors. The geometric
2840  factors, G_OTHER, varies with atmospheric layer number for
2841  aircraft observational geometries.
2842  */
2843  for (i = 0; i < nl; i++)
2844  vmrm[i] = sclgas * (float) 355.0 * (float) 1.e-6 * geometry3_.g_other[i];
2845 
2846  apply_gas_trans(k_surf, levels, bands, abscf_co2, dp, vmrm, init_speccal1_.tran_hi_others);
2847  // for (i = 0; i < bands; i++)
2848  // if (init_speccal1_.tran_hi_others[i]<1.0) printf("CO2: %d tran_hi_others=%f\n",i, init_speccal1_.tran_hi_others[i]);
2849 
2850  }
2851 
2852 
2853  /*
2854  --------------------------------------------
2855  For N2O transmittance calculation.
2856  */
2857  if (n2o) {
2858  sclgas = 0.3;
2859  for (i = 0; i < nl; i++)
2860  vmrm[i] = sclgas * (float) 1.e-6 * geometry3_.g_other[i];
2861 
2862  apply_gas_trans(k_surf, levels, bands, abscf_n2o, dp, vmrm, init_speccal1_.tran_hi_others);
2863  // for (i = 0; i < bands; i++)
2864  // if (init_speccal1_.tran_hi_others[i]<1.0)printf("N2O: %d tran_hi_others=%f\n",i, init_speccal1_.tran_hi_others[i]);
2865  }
2866  /*
2867  --------------------------------------------
2868  For CO transmittance calculation.
2869  */
2870  if (co) {
2871  sclgas = 0.1;
2872  for (i = 0; i < nl; i++)
2873  vmrm[i] = sclgas * (float) 1.e-6 * geometry3_.g_other[i];
2874 
2875  apply_gas_trans(k_surf, levels, bands, abscf_co, dp, vmrm, init_speccal1_.tran_hi_others);
2876  // for (i = 0; i < bands; i++)
2877  // if (init_speccal1_.tran_hi_others[i]<1.0)printf("CO: %d tran_hi_others=%f\n",i, init_speccal1_.tran_hi_others[i]);
2878  }
2879  /*
2880 --------------------------------------------
2881  For CH4 transmittance calculation.
2882  For assigning CH4 VMRM
2883  NOTE: The scaling factor of 0.8 for the CH4 VMRS was obtained by comparing
2884  transmittance spectra calculated using our program, which is based on
2885  the Malkmus narrow band spectral model, with a ratioed spectrum
2886  provided by G. C. Toon at Jet Propulsion Laboratory (JPL). The JPL
2887  ratio spectrum was obtained by ratioing a high resolution (0.005 cm-1)
2888  solar spectrum measured at ground level against a solar spectrum
2889  measured at an altitude of approximately 35 km with the same Fourier
2890  Transform spectrometer. The high resolution ratio spectrum was
2891  degraded to a resolution of 10 nm during our derivation of the
2892  scaling factor for the CH4 VMRS.
2893 
2894  sclch4=0.8;
2895  */
2896  if (ch4) {
2897  sclgas = 1.0;
2898  for (i = 0; i < nl; i++)
2899  vmrm[i] = sclgas * (float) 1.6 * (float) 1.0e-6 * geometry3_.g_other[i]; //I don't know why sclgas=1 and then it's multiplied by 1.6 - rjh
2900 
2901  apply_gas_trans(k_surf, levels, bands, abscf_ch4, dp, vmrm, init_speccal1_.tran_hi_others);
2902  // for (i = 0; i < bands; i++)
2903  // if (init_speccal1_.tran_hi_others[i]<1.0)printf("CH4: %d tran_hi_others=%f\n",i, init_speccal1_.tran_hi_others[i]);
2904 
2905  }
2906  if (o2) {
2907  /*
2908  * ***Modified by Bo-Cai Gao on 2/7/2013 to increase O2 absorption
2909  --- coefficients by the factor SCL_O2 for wavelengths > 1.2 micron
2910  --- in order to model properly the atmospheric O2 band centered
2911  --- near 1.265 micron.
2912  *
2913  */
2914  sclgas = 0.21;
2915 
2916  if ((sumcf = (float *) calloc(bands, sizeof (float))) == NULL) {
2917  printf("-E- : Error allocating memory to sumcf\n");
2918  exit(FATAL_ERROR);
2919  }
2920 
2921  for (i = 0; i < nl; i++)
2922  vmrm[i] = sclgas * geometry3_.g_other[i];
2923 
2924  for (i = 0; i < bands; i++)
2925  sumcf[i] = 0.;
2926  /*
2927  * ***Modified by Bo-Cai Gao on 2/7/2013 to increase O2 absorption
2928  --- coefficients by the factor SCL_O2 for wavelengths > 1.2 micron
2929  --- i& < 1.3333 micron in order to model properly the atmospheric
2930  --- O2 band centered near 1.265 micron.
2931  *
2932  */
2933  sclgas = 2.60;
2934  //printf("ATREM: INIT_SPECCAL - be sure to change below line to i=k_surf and i - k_surf (NOT i - k_surf + 1) once final conversion to C is complete\n ");
2935 
2936  for (i = k_surf - 1; i < levels; i++)
2937  for (j = 0; j < bands; j++)
2938  sumcf[j] -= abscf_o2[bands * i + j] * dp[i - k_surf + 1] * vmrm[i - k_surf + 1];
2939 
2940  for (i = 0; i < 9000; i++)
2941  init_speccal1_.tran_hi_others[i] *= exp(sumcf[i] * Q * (float) 28.966 / (float) 6.0225E+23 / (float) 1.0E-06);
2942  for (i = 9000; i < 106600; i++)
2943  init_speccal1_.tran_hi_others[i] *= exp(sumcf[i] * Q * sclgas * (float) 28.966 / (float) 6.0225E+23 / (float) 1.0E-06);
2944  for (i = 106600; i < NP_HI; i++)
2945  init_speccal1_.tran_hi_others[i] *= exp(sumcf[i] * Q * (float) 28.966 / (float) 6.0225E+23 / (float) 1.0E-06);
2946  // for (i = 0; i < bands; i++) {
2947  // printf("OX: %d tran_hi_others %f\n",i, init_speccal1_.tran_hi_others[i]);
2948  // if (i==94923) {
2949  // printf("OX: %d %f %f\n",i,exp(sumcf[i]*Q*sclgas* (float)28.966 / (float)6.0225E+23 /(float) 1.0E-06),sumcf[i]);
2950  // for (k=k_surf-1; k<levels;k++)
2951  // printf("OX: ABSCF: %d %f %f %f\n",k,abscf_o2[bands * k + i],dp[k-k_surf+1],vmrm[k-k_surf+1]); }
2952  // }
2953  free(sumcf);
2954  }
2955 }
2956 
2957 /*
2958  * Converted to C 4/2016 - R. Healy
2959  * Note this is the merger of two subroutines (TRANCAL with TRAN_TABLE) in the original Atrem Fortran code.
2960  * The original code looped through the 60 Water Vapor values in TRAN_TABLE and called TRANCAL.
2961 !********************************************************************************
2962 !* *
2963 !* Name: TRAN_TABLE *
2964 !* Purpose: This subroutine generates a table consisting of 60 atmospheric *
2965 !* transmittance spectra at the solar and observational *
2966 !* geometry and with 60 column water vapor values. The table also *
2967 !* includes the total amounts of column water vapor used in the *
2968 !* calculations, and the 3-channel ratios calculated from the window *
2969 !* and absorption channels in and around the 0.94- and 1.14-um water *
2970 !* vapor bands. *
2971 !* Parameters: none *
2972 !* Algorithm: For each of the 60 water vapor amounts, calculate the *
2973 !* atmospheric transmittance, and save *
2974 !* Globals Used: NH2O - number of column water vapor values *
2975 !* VAPTT - geometrically adjusted water vapor total *
2976 !* R094 - channel ratio for .94 um region *
2977 !* R114 - channel ratio for 1.14 um region *
2978 !* TRNCAL - atmospheric transmittance spectra *
2979 !* Global Output: VAPTOT() - array containing geometrically adjusted water *
2980 !* vapor values *
2981 !* ROP94() - array containing channel ratios for .94 um region *
2982 !* R1P14() - array containing channel ratios for 1.14 um region*
2983 !* TRNTBL() - 2 dimensional array containing one transmittance *
2984 !* spectrum for each column water vapor amount *
2985 !* Return Codes: none *
2986 !* Special Considerations: none *
2987 !* *
2988 !********************************************************************************
2989 !********************************************************************************
2990 !* TRANCAL combined with TRAN_TABLE by R. Healy 4/28/2015 *
2991 !* Name: TRANCAL *
2992 !* Purpose: This program calculates combined transmittances of H2O, CO2, O3, *
2993 !* N2O, CO, CH4, and O2. *
2994 !* Parameters: none. *
2995 !* Algorithm: The calculations were based on the line-by-line absorption *
2996 !* parameters supplied by William R. Ridgway of NASA/GSFC. *
2997 !* Global output:VAPTT - geometrically adjusted water vapor total. *
2998 !* R094 - channel ratio for 0.94 um region. *
2999 !* R114 - channel ratio for 1.14 um region. *
3000 !* TRNCAL - total transmittances of all gases that match the *
3001 !* resolutions of imaging spectrometers. *
3002 !* Return Codes: none. *
3003 !* Special Considerations: The high resolution (0.05 cm-1) line-by-line *
3004 !* absorption parameters cover the 0.555 - 3.33 micron spectral range *
3005 !* (3000 - 18000 cm-1). The medium resolution ozone absorption *
3006 !* coefficients covers the 0.3-0.8 um spectral range. The line-by-line *
3007 !* high resolution spectra were first smoothed to medium resolution *
3008 !* spectra (resolution = 0.2 nm, wavelength spacing = 0.1 nm) covering *
3009 !* the 0.56 - 3.1 micron spectral region. The medium resolution spectra *
3010 !* of O3 and other gases are combined (in SUBROUTINE TRAN_SMOOTH) to form *
3011 !* a single medium resolution spectrum from 0.3 to 3.1 micron. This *
3012 !* combined spectrum (medium resolution) is then smoothed to lower *
3013 !* resolutions to match the resolutions of imaging spectrometers. The *
3014 !* smoothing is also done in SUBROUTINE TRAN_SMOOTH. *
3015 !* *
3016 !********************************************************************************
3017  */
3018 
3019 void tran_table() {
3020 
3021  int32_t co2 = getinput1_.co2, o2 = getinput1_.o2, n2o = getinput1_.n2o, co = getinput1_.co, ch4 = getinput1_.ch4, no2 = getinput1_.no2, o3 = getinput1_.o3;
3022  float *wavobs = getinput4_.wavobs, *dp = init_speccal5_.dp, *vmrm = init_speccal5_.vmrm;
3023  float *vmr = getinput3_.vmr;
3024  int32_t nbands = getinput5_.nbands, nh2o = init_speccal3_.nh2o, fullcalc = getinput5_.full_calc;
3025  int32_t nl = getinput3_.nl, splitpaths = geometry_l2gen_.splitpaths, ja = geometry_l2gen_.ja, jb = geometry_l2gen_.jb;
3026  float Q = model_adj1_.q;
3027  float mu = geometry5_.mu, mu0 = geometry5_.mu0, *ssh2o = geometry2_.ssh2o, *g_vap = geometry3_.g_vap;
3028  float dvap_plane = model_adj3_.dvap_plane, dvap_layer = model_adj3_.dvap_layer;
3029  int32_t k_plane = model_adj3_.k_plane, k_surf = model_adj4_.k_surf;
3030  float *vaptot = tran_table1_.vaptot;
3031  int i, j, k, n;
3032  static int firstCall = 1;
3033 
3034  static float *abscf_h2o;
3035  static float *tran_hi, *wavno_hi;
3036 
3037  char *filedir;
3038  char *abscf_file = "abscf_gas.nc";
3039  char filename[FILENAME_MAX];
3040  size_t start[3], count[3];
3041  int xid, yid;
3042  int retval, sds_id;
3043  static size_t bands, levels;
3044 
3045  static int32_t n_g = 7, doSmooth;
3046  float tg[7] = {0.0, 0.379, 0.5106, 0.81, 0.9548, 0.9933, 1.0}, f1, f2;
3047 
3048  static float *tkcdf;
3049  static float *sum_kd, *sumcf, **sumcf_s;
3050  float vmrm_s[2][MODELMAX];
3051 
3052  if (firstCall) {
3053 
3054  if (fullcalc)
3055  printf("ATREM: **WARNING** Full_calc is on. Processing will be extremely slow...\n");
3056  else
3057  printf("ATREM: Full_calc is off. Processing using k-distribution method...\n");
3058 
3059  doSmooth = co2 || n2o || co || ch4 || o2 || o3 || no2;
3060 
3061  // Open the netcdf4 input file
3062  if ((filedir = getenv("OCDATAROOT")) == NULL) {
3063  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
3064  exit(FATAL_ERROR);
3065  }
3066 
3067  strcpy(filename, filedir);
3068  strcat(filename, "/common/");
3069  // strcpy(filename,input_l2gen_.filename);
3070  strcat(filename, abscf_file);
3071  strcat(filename, "\0");
3072  retval = nc_open(filename, NC_NOWRITE, &sds_id);
3073  if (retval != NC_NOERR) {
3074  fprintf(stderr, "-E- %s line %d: nc_open(%s) failed.\n",
3075  __FILE__, __LINE__, filename);
3076  exit(FATAL_ERROR);
3077  }
3078  // Now read in the absorption coefficients
3079  retval = nc_inq_dimid(sds_id, "level", &xid)
3080  || nc_inq_dimlen(sds_id, xid, &levels)
3081  || nc_inq_dimid(sds_id, "band", &yid)
3082  || nc_inq_dimlen(sds_id, yid, &bands);
3083 
3084  if (retval) {
3085  fprintf(stderr, "-E- %s line %d: nc_inq_dimid(%s) failed. %d %d %d \n",
3086  __FILE__, __LINE__, filename, (int) levels, (int) bands, retval);
3087  exit(FATAL_ERROR);
3088  }
3089 
3090  if (!(abscf_h2o = (float *) malloc(levels * bands * sizeof (float)))) {
3091  printf("-E- : Error allocating memory to abscf_h2o\n");
3092  exit(FATAL_ERROR);
3093  }
3094 
3095 
3096  start[0] = 0;
3097  start[1] = 0;
3098  start[2] = 0;
3099  count[0] = levels;
3100  count[1] = bands;
3101  count[2] = 0;
3102 
3103  retval = nc_inq_varid(sds_id, "abscf_h2o", &xid);
3104 
3105  if (retval != NC_NOERR) {
3106  fprintf(stderr,
3107  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
3108  __FILE__, __LINE__, filename, "abscf_h2o");
3109  exit(FATAL_ERROR);
3110  }
3111  retval = nc_get_vara_float(sds_id, xid, start, count, abscf_h2o);
3112  if (retval != NC_NOERR) {
3113  fprintf(stderr,
3114  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
3115  __FILE__, __LINE__, filename, "abscf_h2o");
3116  exit(FATAL_ERROR);
3117  }
3118 
3119  if (!wavno_hi && !(wavno_hi = (float *) malloc(NP_HI * sizeof (float *)))) {
3120  printf("-E- : Error allocating memory to wavno_hi\n");
3121  exit(FATAL_ERROR);
3122  }
3123 
3124  retval = nc_inq_varid(sds_id, "waveno", &xid);
3125 
3126  count[0] = bands;
3127  count[1] = 0;
3128 
3129  retval = nc_get_vara_float(sds_id, xid, start, count, wavno_hi);
3130  if (retval != NC_NOERR) {
3131  fprintf(stderr,
3132  "-E- %s line %d: nc_get_vara_float failed for file, %s field, %s.\n",
3133  __FILE__, __LINE__, filename, "wavno_hi");
3134  exit(FATAL_ERROR);
3135  }
3136 
3137  nc_close(sds_id);
3138 
3139  if (!(tkcdf = (float *) calloc(n_g * nbands * levels, sizeof (float)))) {
3140  printf("-E- : Error allocating memory to tkcdf\n");
3141  exit(FATAL_ERROR);
3142  }
3143 
3144  if (!(sum_kd = (float *) calloc(nbands, sizeof (float)))) {
3145  printf("-E- : Error allocating memory to sum_kd\n");
3146  exit(FATAL_ERROR);
3147  }
3148 
3149  }
3150 
3151  /*
3152  ! Moved from INIT_SPECAL because it seems more appropriate here. - rjh 01/04/2016
3153  !C Initial water vapor VMRs for repeated use in other subroutines and
3154  !C adjust layered water vapor VMRM with geometrical factors.
3155  */
3156  //if (firstCall) printf("ATREM: TRAN_CAL - be sure to remove below 2 lines (k_surf--;k_plane--) once final conversion to C is complete\n ");
3157  k_plane--; // RJH - this is a relic of the fortran code - may want change this throughout the code later
3158  k_surf--; // Ditto
3159 
3160  for (i = 0; i < nl; i++) {
3161  vmrm[i] = (vmr[i] + vmr[i + 1]) / (float) 2.0;
3162  /*
3163  Scale the VMRM by the two-way path geometrical factors. The geometric
3164  factors, G_VAP, varies with atmospheric layer number for
3165  aircraft observational geometries.
3166  */
3167  if (splitpaths) { // perform the split path calculations
3168  vmrm_s[0][i] = vmrm[i] * mu0;
3169  if (i < k_plane)
3170  vmrm_s[1][i] = vmrm[i] * mu;
3171  else {
3172  if (i == k_plane)
3173  vmrm_s[1][i] = vmrm[i] * mu * dvap_plane / dvap_layer;
3174  else
3175  vmrm_s[1][i] = 0.;
3176 
3177  }
3178  }
3179  vmrm[i] *= g_vap[i];
3180  }
3181 
3182  if (firstCall || fullcalc) {
3183  // For each water vapor amount, calculate the geometrically adjusted water
3184  // vapor amount, the channel ratios, and the transmittance spectrum.
3185  if (!sumcf && !(sumcf = (float *) calloc(NP_HI, sizeof (float)))) { //FREE
3186  printf("-E- : Error allocating memory to sumcf\n");
3187  exit(FATAL_ERROR);
3188  }
3189  if (splitpaths) {
3190  if (!sumcf_s) {
3191  if (!(sumcf_s = (float **) malloc(2 * sizeof (float*)))) { //FREE
3192  printf("-E- : Error allocating memory to sumcf_s\n");
3193  exit(FATAL_ERROR);
3194  }
3195  for (i = 0; i < 2; i++) {
3196  if (!(sumcf_s[i] = (float *) calloc(NP_HI, sizeof (float)))) { //FREE
3197  printf("-E- : Error allocating memory to sumcf_s\n");
3198  exit(FATAL_ERROR);
3199  }
3200  }
3201  }
3202  }
3203  for (i = 0; i < bands; i++) {
3204  sumcf[i] = 0;
3205  if (splitpaths) {
3206  sumcf_s[0][i] = 0;
3207  sumcf_s[1][i] = 0;
3208  }
3209  for (j = k_surf; j < levels; j++) {
3210  sumcf[i] -= abscf_h2o[bands * j + i] * dp[j - k_surf] * vmrm[j - k_surf];
3211  // printf("sumcf[%d]=%g abscf_h2o[%d][%d]=%g %f %f \n ",i,sumcf[i],j,i,abscf_h2o[bands*j+i],dp[j-k_surf],vmrm[j-k_surf]);
3212  if (splitpaths) {
3213  sumcf_s[0][i] -= abscf_h2o[bands * j + i] * dp[j - k_surf] * vmrm_s[0][j - k_surf];
3214  sumcf_s[1][i] -= abscf_h2o[bands * j + i] * dp[j - k_surf] * vmrm_s[1][j - k_surf];
3215  }
3216  }
3217  }
3218 
3219  if (!tran_hi && !(tran_hi = (float *) malloc(nh2o * NP_HI * sizeof (float *)))) {
3220  printf("-E- : Error allocating memory to tran_hi\n");
3221  exit(FATAL_ERROR);
3222  }
3223 
3224  for (i = 0; i < nh2o; i++) {
3225  for (j = 0; j < bands; j++) {
3226  tran_hi[i * NP_HI + j] = exp(sumcf[j] * ssh2o[i] * Q * (float) 28.966 / (float) 6.0225e+23 / (float) 1.0e-6);
3227  }
3228 
3229  // Total amount of water vapor (in unit of cm) corresponding to the spectrum.
3230 
3231  vaptot[i] = geometry4_.vap_slant_mdl * ssh2o[i];
3232  }
3233 
3234  if (splitpaths) {
3235  for (j = 0; j < 2; j++) {
3236  for (i = 0; i < bands; i++) {
3237  // printf("TRAN_TABLE: SSH2O_S JA: %d %g \n", ja, geometry5_.ssh2o_s[j][ja]);
3238  // printf("TRAN_TABLE: SSH2O_S JA+1: %d %g \n",ja+1,geometry5_.ssh2o_s[j][ja+1]);
3239  // printf("TRAN_TABLE: SSH2O_S JB: %d %g \n", jb, geometry5_.ssh2o_s[j][jb]);
3240  // printf("TRAN_TABLE: SSH2O_S JB+1: %d %g \n",jb+1,geometry5_.ssh2o_s[j][jb+1]);
3241  tran_tables_.tran_hi_sa [j][i] = exp(sumcf_s[j][i] * geometry5_.ssh2o_s[j][ja - 1] * Q * (float) 28.966 / (float) 6.0225e+23 / (float) 1.0e-6);
3242  tran_tables_.tran_hi_sap1[j][i] = exp(sumcf_s[j][i] * geometry5_.ssh2o_s[j][ja ] * Q * (float) 28.966 / (float) 6.0225e+23 / (float) 1.0e-6);
3243  tran_tables_.tran_hi_sb [j][i] = exp(sumcf_s[j][i] * geometry5_.ssh2o_s[j][jb - 1] * Q * (float) 28.966 / (float) 6.0225e+23 / (float) 1.0e-6);
3244  tran_tables_.tran_hi_sbp1[j][i] = exp(sumcf_s[j][i] * geometry5_.ssh2o_s[j][jb ] * Q * (float) 28.966 / (float) 6.0225e+23 / (float) 1.0e-6);
3245  }
3246 
3247  }
3248  }
3249  }
3250 
3251  if (!fullcalc) {
3252  // Calculate the k-distrubition coefficients
3253 
3254  if (firstCall) {
3255  //kdist_gas_abs_(tkcdf,abscf_h2o,&nphi,wavno_hi,wavobs,&nbands) ; //FORTRAN
3256  kdistgasabs(tkcdf, abscf_h2o, wavno_hi, wavobs, bands, levels, nbands); //C
3257  }
3258  // Perform the k-distribution calculation Trapezoidal integral for transmittance over the NBANDS bands
3259 
3260  for (i = 0; i < nh2o; i++) {
3261  for (j = 0; j < nbands; j++) tran_table1_.tran_kd[i][j] = 1.0;
3262 
3263  for (j = k_surf; j < levels; j++) {
3264 
3265  //tkcdf[n_g][nbands][levels]
3266  for (n = 0; n < nbands; n++) {
3267  sum_kd[n] = 0.0;
3268  f1 = exp(-tkcdf[0 * nbands * levels + levels * n + j] * dp[j - k_surf] * vmrm[j - k_surf] * ssh2o[i]);
3269 
3270  for (k = 0; k < n_g - 1; k++) {
3271  f2 = exp(-tkcdf[(k + 1) * nbands * levels + levels * n + j] * dp[j - k_surf] * vmrm[j - k_surf] * ssh2o[i]);
3272  sum_kd[n] += (f1 + f2)*(tg[k + 1] - tg[k]) / (float) 2.0;
3273  f1 = f2;
3274  }
3275  tran_table1_.tran_kd[i][n] *= sum_kd[n];
3276  }
3277  }
3278  // Total amount of water vapor (in unit of cm) corresponding to the spectrum.
3279 
3280  vaptot[i] = geometry4_.vap_slant_mdl * ssh2o[i];
3281  }
3282  }
3283 
3284  if (fullcalc || firstCall) {
3285 
3286  // Smooth the high resolution spectra to resolutions of measured spectrum
3287  //tran_smooth_(tran_hi); //FORTRAN
3288  tran_smooth(tran_hi);
3289  if (!fullcalc) {
3290  free(wavno_hi);
3291  free(tran_hi);
3292  }
3293 
3294  }
3295 
3296  for (i = 0; i < nbands; i++) {
3297  tran_table1_.trntblo[i] = 1.0;
3298  tran_table_l2gen_.tg_solo[i] = 1.0;
3299  tran_table_l2gen_.tg_seno[i] = 1.0;
3300  }
3301 
3302  // if (doSmooth) tran_smooth_others_(); //FORTRAN
3303  if (doSmooth) tran_smooth_others();
3304 
3305  if (!fullcalc) {
3306 
3307  if (firstCall)
3308  for (k = 0; k < nbands; k++)
3309  printf("%d ATREM: FAST)TRNTBL=%f %f\n", k, tran_table1_.tran_kd[59][k] * tran_table1_.diff_tran[0][k], tran_table1_.trntblo[k]);
3310 
3311  for (i = 0; i < nh2o; i++) {
3312  for (k = 0; k < nbands; k++) {
3313  /*
3314  Multiplying the sensor resolution water vapor transmittance with combined
3315  transmittances of CO2, N2O, CO, CH4, and O2:
3316  */
3317  tran_table1_.trntbl[i][k] = tran_table1_.tran_kd[i][k] * tran_table1_.diff_tran[i][k] * tran_table1_.trntblo[k];
3318  if (firstCall)
3319  printf("ATREM: DIFF_TRAN: %d %d %f %f %f %f\n", i, k,
3320  tran_table1_.trntbl[i][k], tran_table1_.trntblo[k], tran_table1_.tran_kd[i][k], tran_table1_.diff_tran[i][k]);
3321  }
3322  }
3323 
3324  } else {
3325 
3326  /*
3327  Multiplying the sensor resolution water vapor transmittance with combined
3328  transmittances of CO2, N2O, CO, CH4, and O2:
3329  */
3330  for (k = 0; k < nbands; k++) {
3331  for (i = 0; i < nh2o; i++) {
3332  tran_table1_.trntbl[i][k] *= tran_table1_.trntblo[k];
3333  }
3334 
3335  if (splitpaths) {
3336  tran_table_l2gen_.tg_sol[k] *= tran_table_l2gen_.tg_solo[k];
3337  tran_table_l2gen_.tg_sen[k] *= tran_table_l2gen_.tg_seno[k];
3338  }
3339 
3340  }
3341  }
3342 
3343  //chnlratio_(); // FORTRAN
3344  channelRatio(); // C version
3345 
3346  firstCall = 0;
3347 
3348 
3349 
3350 }
3351 
3352 /* Converted to C from Fortran by R. Healy (richard.healy@nasa.gov) 4/2016
3353  *
3354 !********************************************************************************
3355 !* *
3356 !* Name: TRAN_SMOOTH *
3357 !* Purpose: This program is to smooth the line-by-line high resolution *
3358 !* spectrum to lower resolution spectrum that matches the resolutions *
3359 !* of imaging spectrometer data. *
3360 !* Parameters: none. *
3361 !* Algorithm: The smoothing is done in two stages. The 1st stage is to smooth *
3362 !* the high resolution spectrum to medium resolution spectrum at a *
3363 !* constant FWHM (0.2 nm) and a constant wavelength interval *
3364 !* (0.1 nm). The 2nd stage smoothing is to smooth the medium *
3365 !* resolution spectrum to resolutions of input imaging spectrometer *
3366 !* data. *
3367 !* Globals used: The global variables used are contained in the file *
3368 !* "COMMONS_INC" *
3369 !* Global output: TRNCAL - total transmittances of all gases that match the *
3370 !* resolutions of imaging spectrometers. *
3371 !* Return Codes: none. *
3372 !* *
3373 !********************************************************************************
3374  */
3375 void tran_smooth(float *tran_hi) {
3376 
3377  static int firstCall = 1;
3378 
3379  /*
3380  * Converted to C from Fortran by R. Healy (richard.healy@nasa.gov) 4/2016
3381  *
3382  First stage of smoothing - smooth line-by-line high resolution spectrum with
3383  over 300,000 points (point spacing of 0.05 cm-1) to medium resolution
3384  spectrum (resolution of 0.2 nm and point spacing of 0.1 nm) with about
3385  25,000 points.
3386 
3387  The smoothing of line-by-line spectrum is done in wavenumber domain. For
3388  a spectrum with a constant 0.2 nm resolution in wavelength domain, it has
3389  variable resolution in wavenumber domain. This effect is properly taken
3390  care of in the design of smoothing functions.
3391 
3392  Because the high resolution spectrum is in wavenumber units (cm-1), while
3393  the medium resolution spectrum is in wavelength units, the two kinds of
3394  grids do not automatically match. In order to match the grids, arrays
3395  of INDEX_MED and TRAN_MED_INDEX are specially designed. The desired
3396  medium resolution spectrum, TRAN_MED, at constant 0.1 nm point spacing
3397  and 0.2 nm resolution is obtained through linear interpolation of
3398  TRAN_MED_INDEX array.
3399 
3400  */
3401  int32_t i, j, k;
3402  float *wavln_med_index = init_speccal13_.wavln_med_index, *fwhm = getinput4_.fwhm, *wavln_std = init_speccal12_.wavln_std;
3403  float *wavobs = getinput4_.wavobs;
3404  static float **finstr_wavno, *fwhm_wavno;
3405  static int32_t ncvhf_wavno[NP_MED], *ia, ia_p1;
3406  int32_t ncvtot_wavno, ncvtot, nbands = getinput5_.nbands, *ncvhf = init_speccal15_.ncvhf, *index_med = init_speccal13_.index_med;
3407 
3408  float *tran_med_sa_sol = tran_tables1_.tran_med_sa_sol, *tran_med_sap1_sol = tran_tables1_.tran_med_sap1_sol;
3409  float *tran_med_sa_sen = tran_tables1_.tran_med_sa_sen, *tran_med_sap1_sen = tran_tables1_.tran_med_sap1_sen;
3410  float *tran_med_sb_sol = tran_tables1_.tran_med_sb_sol, *tran_med_sbp1_sol = tran_tables1_.tran_med_sbp1_sol;
3411  float *tran_med_sb_sen = tran_tables1_.tran_med_sb_sen, *tran_med_sbp1_sen = tran_tables1_.tran_med_sbp1_sen;
3412  float *tran_med_index_sa_sol = tran_tables1_.tran_med_index_sa_sol, *tran_med_index_sap1_sol = tran_tables1_.tran_med_index_sap1_sol;
3413  float *tran_med_index_sa_sen = tran_tables1_.tran_med_index_sa_sen, *tran_med_index_sap1_sen = tran_tables1_.tran_med_index_sap1_sen;
3414  float *tran_med_index_sb_sol = tran_tables1_.tran_med_index_sb_sol, *tran_med_index_sbp1_sol = tran_tables1_.tran_med_index_sbp1_sol;
3415  float *tran_med_index_sb_sen = tran_tables1_.tran_med_index_sb_sen, *tran_med_index_sbp1_sen = tran_tables1_.tran_med_index_sbp1_sen;
3416  float *tran_std_sa_sol = tran_tables1_.tran_std_sa_sol, *tran_std_sap1_sol = tran_tables1_.tran_std_sap1_sol;
3417  float *tran_std_sa_sen = tran_tables1_.tran_std_sa_sen, *tran_std_sap1_sen = tran_tables1_.tran_std_sap1_sen;
3418  float *tran_std_sb_sol = tran_tables1_.tran_std_sb_sol, *tran_std_sbp1_sol = tran_tables1_.tran_std_sbp1_sol;
3419  float *tran_std_sb_sen = tran_tables1_.tran_std_sb_sen, *tran_std_sbp1_sen = tran_tables1_.tran_std_sbp1_sen;
3420  float *tg_sol = tran_table_l2gen_.tg_sol, *tg_sen = tran_table_l2gen_.tg_sen;
3421 
3422  float tran_ia[NH2OMAX], tran_iap1[NH2OMAX];
3423  float tran_ia_sa_sol, tran_ia_sa_sen, tran_ia_sap1_sol, tran_ia_sap1_sen;
3424  float tran_iap1_sa_sol, tran_iap1_sa_sen, tran_iap1_sap1_sol, tran_iap1_sap1_sen;
3425  float tran_ia_sb_sol, tran_ia_sb_sen, tran_ia_sbp1_sol, tran_ia_sbp1_sen;
3426  float tran_iap1_sb_sol, tran_iap1_sb_sen, tran_iap1_sbp1_sol, tran_iap1_sbp1_sen;
3427  float sumins, dlt, fj, fjm1, dlt_ia, fia, fia_p1;
3428  float **tran_std, **tran_med;
3429  float tg_sol_a, tg_sen_a, tg_sol_ap1, tg_sen_ap1;
3430  float tg_sol_b, tg_sen_b, tg_sol_bp1, tg_sen_bp1;
3431  float f1a = geometry_l2gen_.f1a, f2a = geometry_l2gen_.f2a;
3432  float f1b = geometry_l2gen_.f1b, f2b = geometry_l2gen_.f2b;
3433  int splitpaths = geometry_l2gen_.splitpaths;
3434  int32_t np_std = NP_STD, ndx1, ndx2;
3435 
3436  if (firstCall) {
3437  if (!(fwhm_wavno = (float *) malloc(NP_MED * sizeof (float)))) {
3438  printf("-E- : Error allocating memory to fwhm_wavno\n");
3439  exit(FATAL_ERROR);
3440  }
3441  if (!(finstr_wavno = (float **) malloc(NP_MED * sizeof (float*)))) {
3442  printf("-E- : Error allocating memory to finstr_wavno\n");
3443  exit(FATAL_ERROR);
3444  }
3445  for (j = 0; j < NP_MED; j++) {
3446  if (!(finstr_wavno[j] = (float *) malloc(NINSTRF * sizeof (float)))) {
3447  printf("-E- : Error allocating memory to finstr_wavno\n");
3448  exit(FATAL_ERROR);
3449  }
3450  fwhm_wavno[j] = (float) 10000. * DLT_MED / (wavln_med_index[j] * wavln_med_index[j]);
3451  ncvhf_wavno[j] = (FACDLT * fwhm_wavno[j] / DWAVNO + 1.);
3452  ncvtot_wavno = 2 * ncvhf_wavno[j] - 1;
3453 
3454  sumins = 0;
3455  for (i = ncvhf_wavno[j] - 1; i < ncvtot_wavno; i++) {
3456  finstr_wavno[j][i] = exp(-CONST1 * pow((float) (i - ncvhf_wavno[j] + 1) * DWAVNO / fwhm_wavno[j], 2));
3457  sumins += finstr_wavno[j][i];
3458  // if (j==0) printf("A:finstr_wavno sumins: %d %g %g %d\n",i,sumins,finstr_wavno[j][i],i-ncvhf_wavno[j]+1);
3459  }
3460 
3461  for (i = 0; i < ncvhf_wavno[j] - 1; i++) {
3462  finstr_wavno[j][i] = finstr_wavno[j][ncvtot_wavno - i - 1];
3463  sumins += finstr_wavno[j][i];
3464  // if (j==0) printf("B:finstr_wavno sumins: %d %g %g %d\n",i,sumins,finstr_wavno[j][i],ncvtot_wavno-i-1);
3465  }
3466 
3467  sumins *= DWAVNO;
3468 
3469  for (i = 0; i < ncvtot_wavno; i++) {
3470  finstr_wavno[j][i] *= (DWAVNO / sumins);
3471  // printf("RJH: finstr_wavno %d %d %g %g \n",i,j,finstr_wavno[j][i],sumins);
3472  }
3473 
3474  }
3475 
3476  ia = (int32_t*) allocateMemory(nbands * sizeof (int32_t), "ia");
3477 
3478  // Index searching...
3479  //
3480  //Calculate instrumental response functions...
3481 
3482 
3483  for (j = 0; j < nbands; j++) {
3484  sumins = 0.0;
3485  ncvtot = 2 * ncvhf[j] - 1;
3486  for (i = ncvhf[j] - 1; i < ncvtot; i++) {
3487  init_speccal15_.finstr[j][i] = exp(-CONST1 * pow((float) (i - ncvhf[j] + 1) * DWAVLN / fwhm[j], 2));
3488  sumins += init_speccal15_.finstr[j][i];
3489  // if (j==0) printf("c:A:finstr sumins: %d %g %g %d fwhm=%g %g %g %g %g %g\n",i,sumins,init_speccal15_.finstr[j][i],i-ncvhf[j]+1,fwhm[j],DWAVLN,
3490  // (float) (i-ncvhf[j]+1)*DWAVLN/fwhm[j],pow( (float) (i-ncvhf[j]+1)*DWAVLN/fwhm[j],2),-CONST1*pow( (float) (i-ncvhf[j]+1)*DWAVLN/fwhm[j],2),
3491  // exp(-CONST1*pow( (float) (i-ncvhf[j]+1)*DWAVLN/fwhm[j],2)));
3492  }
3493 
3494  for (i = 0; i < ncvhf[j] - 1; i++) {
3495  init_speccal15_.finstr[j][i] = init_speccal15_.finstr[j][ncvtot - i - 1];
3496  sumins += init_speccal15_.finstr[j][i];
3497  // if (j==0) printf("c:B:finstr sumins: %d %g %g %d\n",i,sumins,init_speccal15_.finstr[j][i],ncvtot-i-1);
3498  }
3499 
3500  sumins *= DWAVLN;
3501 
3502  for (i = 0; i < ncvtot; i++)
3503  init_speccal15_.finstr[j][i] *= (DWAVLN / sumins);
3504 
3505  ia[j] = hunt(wavln_std, np_std, (double) wavobs[j], ia[j]);
3506  }
3507 
3508  firstCall = 0;
3509  }
3510 
3511  if (!(tran_std = (float **) malloc(NH2OMAX * sizeof (float *)))) {
3512  printf("-E- : Error allocating memory to tran_std\n");
3513  exit(FATAL_ERROR);
3514  }
3515  if (!(tran_med = (float **) malloc(NH2OMAX * sizeof (float *)))) {
3516  printf("-E- : Error allocating memory to tran_med\n");
3517  exit(FATAL_ERROR);
3518  }
3519 
3520  for (i = 0; i < NH2OMAX; i++) {
3521  if (!(tran_std[i] = (float *) calloc(NP_STD, sizeof (float *)))) {
3522  printf("-E- : Error allocating memory to tran_std\n");
3523  exit(FATAL_ERROR);
3524  }
3525  if (!(tran_med[i] = (float *) calloc(NP_MED, sizeof (float *)))) {
3526  printf("-E- : Error allocating memory to tran_med\n");
3527  exit(FATAL_ERROR);
3528  }
3529  }
3530  for (j = 0; j < NP_MED; j++) {
3531  for (i = 0; i < NH2OMAX; i++)
3532  init_speccal13_.tran_med_index[i][j] = 0.;
3533  if (splitpaths) {
3534  tran_med_index_sa_sen [j] = 0;
3535  tran_med_index_sap1_sen[j] = 0;
3536  tran_med_index_sb_sen [j] = 0;
3537  tran_med_index_sbp1_sen[j] = 0;
3538  tran_med_index_sa_sol [j] = 0;
3539  tran_med_index_sap1_sol[j] = 0;
3540  tran_med_index_sb_sol [j] = 0;
3541  tran_med_index_sbp1_sol[j] = 0;
3542  }
3543  ndx1 = index_med[j]-(ncvhf_wavno[j] - 1);
3544  ndx2 = index_med[j] + ncvhf_wavno[j] - 1;
3545 
3546  for (k = ndx1 - 1; k < ndx2; k++) {
3547  for (i = 0; i < NH2OMAX; i++)
3548  init_speccal13_.tran_med_index[i][j] += tran_hi[i * NP_HI + k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3549 
3550 
3551  // if (i==29)
3552  // printf(" TRAN_MED_INDEX: %d %d <%g> %g %g %d %d %d\n",j,k,init_speccal13_.tran_med_index[i][j],tran_hi[i*NP_HI+k],
3553  // finstr_wavno[j][k-index_med[j]+ncvhf_wavno[j]],
3554  // index_med[j],ncvhf_wavno[j],k-index_med[j]+ncvhf_wavno[j]);
3555  if (splitpaths) {
3556  //Note that the tran_hi_* arrays are 0 for solar path and 1 for sensor
3557  tran_med_index_sa_sol [j] += tran_tables_.tran_hi_sa [0][k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3558  tran_med_index_sa_sen [j] += tran_tables_.tran_hi_sa [1][k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3559  tran_med_index_sap1_sol[j] += tran_tables_.tran_hi_sap1[0][k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3560  tran_med_index_sap1_sen[j] += tran_tables_.tran_hi_sap1[1][k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3561  tran_med_index_sb_sol [j] += tran_tables_.tran_hi_sb [0][k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3562  tran_med_index_sb_sen [j] += tran_tables_.tran_hi_sb [1][k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3563  tran_med_index_sbp1_sol[j] += tran_tables_.tran_hi_sbp1[0][k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3564  tran_med_index_sbp1_sen[j] += tran_tables_.tran_hi_sbp1[1][k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
3565  // if (j==8204) printf("SPLIT: TRAN_MED_INDEX: %d %d <%g> %g %g %d %d %d %g\n",j,k,tran_med_index_sa_sol[j],tran_tables_.tran_hi_sa[0][k],
3566  // finstr_wavno[j][k-index_med[j]+ncvhf_wavno[j]],
3567  // index_med[j],ncvhf_wavno[j],k-index_med[j]+ncvhf_wavno[j],
3568  // tran_tables_.tran_hi_sa [0][k]*finstr_wavno[j][k-index_med[j]+ncvhf_wavno[j]]);
3569  }
3570  }
3571  }
3572  /*
3573  *
3574  Linear interpolation to get TRAN_MED from TRAN_MED_INDEX:
3575  (Note that WAVLN_MED_INDEX(J) >= WAVLN_MED(J) )
3576  *
3577  */
3578  for (i = 0; i < NH2OMAX; i++) {
3579  tran_med[i][0] = init_speccal13_.tran_med_index[i][0];
3580  tran_med[i][NP_MED - 1] = init_speccal13_.tran_med_index[i][NP_MED - 1];
3581  }
3582  if (splitpaths) {
3583  tran_med_sa_sol [0] = tran_med_index_sa_sol [0];
3584  tran_med_sa_sol [NP_MED - 1] = tran_med_index_sa_sol [NP_MED - 1];
3585  tran_med_sa_sen [0] = tran_med_index_sa_sen [0];
3586  tran_med_sa_sen [NP_MED - 1] = tran_med_index_sa_sen [NP_MED - 1];
3587 
3588  tran_med_sap1_sol[0] = tran_med_index_sap1_sol[0];
3589  tran_med_sap1_sol[NP_MED - 1] = tran_med_index_sap1_sol[NP_MED - 1];
3590  tran_med_sap1_sen[0] = tran_med_index_sap1_sen[0];
3591  tran_med_sap1_sen[NP_MED - 1] = tran_med_index_sap1_sen[NP_MED - 1];
3592 
3593  tran_med_sb_sol [0] = tran_med_index_sb_sol [0];
3594  tran_med_sb_sol [NP_MED - 1] = tran_med_index_sb_sol [NP_MED - 1];
3595  tran_med_sb_sen [0] = tran_med_index_sb_sen [0];
3596  tran_med_sb_sen [NP_MED - 1] = tran_med_index_sb_sen [NP_MED - 1];
3597 
3598  tran_med_sbp1_sol[0] = tran_med_index_sbp1_sol[0];
3599  tran_med_sbp1_sol[NP_MED - 1] = tran_med_index_sbp1_sol[NP_MED - 1];
3600  tran_med_sbp1_sen[0] = tran_med_index_sbp1_sen[0];
3601  tran_med_sbp1_sen[NP_MED - 1] = tran_med_index_sbp1_sen[NP_MED - 1];
3602  }
3603 
3604  for (j = 1; j < NP_MED - 1; j++) {
3605  if (wavln_med_index[j] <= init_speccal12_.wavln_med[j]) {
3606  for (i = 0; i < NH2OMAX; i++)
3607  tran_med[i][j] = init_speccal13_.tran_med_index[i][j];
3608  // printf(" TRAN_MED_WAVLN: %d %d <%g> <%g> %g %g\n",j,29,tran_med[29][j],init_speccal13_.tran_med_index[29][j],
3609  // wavln_med_index[j],init_speccal12_.wavln_med[j]);
3610  if (splitpaths) {
3611  tran_med_sa_sol [j] = tran_med_index_sa_sol [j];
3612  tran_med_sa_sen [j] = tran_med_index_sa_sen [j];
3613  tran_med_sap1_sol [j] = tran_med_index_sap1_sol [j];
3614  tran_med_sap1_sen [j] = tran_med_index_sap1_sen [j];
3615  tran_med_sb_sol [j] = tran_med_index_sb_sol [j];
3616  tran_med_sb_sen [j] = tran_med_index_sb_sen [j];
3617  tran_med_sbp1_sol [j] = tran_med_index_sbp1_sol [j];
3618  tran_med_sbp1_sen [j] = tran_med_index_sbp1_sen [j];
3619  }
3620  } else {
3621  dlt = wavln_med_index[j] - wavln_med_index[j - 1];
3622  fjm1 = (wavln_med_index[j] - init_speccal12_.wavln_med[j]) / dlt;
3623  fj = (init_speccal12_.wavln_med [j] - wavln_med_index[j - 1]) / dlt;
3624  for (i = 0; i < NH2OMAX; i++) {
3625  tran_med[i][j] = fjm1 * init_speccal13_.tran_med_index[i][j - 1] + fj * init_speccal13_.tran_med_index[i][j];
3626  // if (i==29)
3627  // printf(" TRAN_MED_INTERP: %d %d <%g> <%g> %g %g %g\n",j,i,init_speccal13_.tran_med_index[i][j-1],init_speccal13_.tran_med_index[i][j],
3628  // fjm1,fj,dlt);
3629  }
3630  if (splitpaths) {
3631  tran_med_sa_sol [j] = fjm1 * tran_med_index_sa_sol [j - 1] + fj * tran_med_index_sa_sol [j];
3632  tran_med_sa_sen [j] = fjm1 * tran_med_index_sa_sen [j - 1] + fj * tran_med_index_sa_sen [j];
3633  tran_med_sap1_sol [j] = fjm1 * tran_med_index_sap1_sol[j - 1] + fj * tran_med_index_sap1_sol[j];
3634  tran_med_sap1_sen [j] = fjm1 * tran_med_index_sap1_sen[j - 1] + fj * tran_med_index_sap1_sen[j];
3635  tran_med_sb_sol [j] = fjm1 * tran_med_index_sb_sol [j - 1] + fj * tran_med_index_sb_sol [j];
3636  tran_med_sb_sen [j] = fjm1 * tran_med_index_sb_sen [j - 1] + fj * tran_med_index_sb_sen [j];
3637  tran_med_sbp1_sol [j] = fjm1 * tran_med_index_sbp1_sol[j - 1] + fj * tran_med_index_sbp1_sol[j];
3638  tran_med_sbp1_sen [j] = fjm1 * tran_med_index_sbp1_sen[j - 1] + fj * tran_med_index_sbp1_sen[j];
3639  // printf("SPLIT: TRAN_MED_SA_SOL= %d <%g> <%g %g> %g %g\n",j,tran_med_sa_sol[j],
3640  // fjm1,fj,tran_med_index_sa_sol[j-1],tran_med_index_sa_sol[j]);
3641  }
3642  }
3643  }
3644 
3645  /*
3646  * --- Here multiplying O3 and NO2 spectra and other spectrum at medium resolution:
3647  *
3648  */
3649  for (j = 0; j < NH2OMAX; j++) {
3650  for (i = 0; i < NPSHIF; i++)
3651  tran_std[j][i] = 1.0;
3652  for (i = NPSHIF; i < NP_STD; i++)
3653  tran_std[j][i] = tran_med[j][i - NPSHIF];
3654 
3655  }
3656  if (splitpaths) {
3657  for (i = 0; i < NP_STD; i++) {
3658  tran_std_sa_sol [i] = 1.0;
3659  tran_std_sa_sen [i] = 1.0;
3660  tran_std_sap1_sol[i] = 1.0;
3661  tran_std_sap1_sen[i] = 1.0;
3662  tran_std_sb_sol [i] = 1.0;
3663  tran_std_sb_sen [i] = 1.0;
3664  tran_std_sbp1_sol[i] = 1.0;
3665  tran_std_sbp1_sen[i] = 1.0;
3666  }
3667  for (i = NPSHIF; i < NP_STD; i++) {
3668  tran_std_sa_sol [i] *= tran_med_sa_sol [i - NPSHIF];
3669  tran_std_sa_sen [i] *= tran_med_sa_sen [i - NPSHIF];
3670  tran_std_sap1_sol[i] *= tran_med_sap1_sol[i - NPSHIF];
3671  tran_std_sap1_sen[i] *= tran_med_sap1_sen[i - NPSHIF];
3672  tran_std_sb_sol [i] *= tran_med_sb_sol [i - NPSHIF];
3673  tran_std_sb_sen [i] *= tran_med_sb_sen [i - NPSHIF];
3674  tran_std_sbp1_sol[i] *= tran_med_sbp1_sol[i - NPSHIF];
3675  tran_std_sbp1_sen[i] *= tran_med_sbp1_sen[i - NPSHIF];
3676  // printf("%d TRAN_STD_SA_SOL= %g %g\n",i,tran_std_sa_sol[i],tran_med_sa_sol[i-NPSHIF]);
3677 
3678  }
3679  }
3680  /*
3681  * The 2nd stage of smoothing - smooth the medium resolution spectrum (resolution
3682  of 0.2 nm and point spacing of 0.1 nm) with about 25,000 points to match
3683  the coarser and variable resolution spectrum from imaging spectrometers.
3684 
3685  Initialize some index parameters:
3686 
3687  */
3688 
3689  for (j = 0; j < nbands; j++) {
3690 
3691  for (i = 0; i < NH2OMAX; i++) {
3692  tran_table1_.trntbl[i][j] = 0.0;
3693  tran_ia[i] = 0.0;
3694  tran_iap1[i] = 0.0;
3695  }
3696 
3697  if (splitpaths) {
3698  tg_sol[j] = 0.0;
3699  tg_sen[j] = 0.0;
3700  tran_ia_sa_sol = 0.0;
3701  tran_ia_sa_sen = 0.0;
3702  tran_ia_sap1_sol = 0.0;
3703  tran_ia_sap1_sen = 0.0;
3704  tran_iap1_sa_sol = 0.0;
3705  tran_iap1_sa_sen = 0.0;
3706  tran_iap1_sap1_sol = 0.0;
3707  tran_iap1_sap1_sen = 0.0;
3708  tran_ia_sb_sol = 0.0;
3709  tran_ia_sb_sen = 0.0;
3710  tran_ia_sbp1_sol = 0.0;
3711  tran_ia_sbp1_sen = 0.0;
3712  tran_iap1_sb_sol = 0.0;
3713  tran_iap1_sb_sen = 0.0;
3714  tran_iap1_sbp1_sol = 0.0;
3715  tran_iap1_sbp1_sen = 0.0;
3716 
3717  }
3718 
3719  // Smoothing ...
3720  for (k = ia[j]-(ncvhf[j] - 1); k < ia[j] + ncvhf[j]; k++) {
3721  for (i = 0; i < NH2OMAX; i++) {
3722  tran_ia[i] += tran_std[i][k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3723  // if(i==29) printf("TRAN_IA: %d %d %d %g %g %d %g %d\n",j,k,ia[j],tran_ia[i],tran_std[i][k],k-ia[j]+ncvhf[j]-1,init_speccal15_.finstr[j][k-ia[j]+ncvhf[j]-1],ncvhf[j]);
3724  }
3725  if (splitpaths) {
3726  tran_ia_sa_sol += tran_std_sa_sol [k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3727  tran_ia_sa_sen += tran_std_sa_sen [k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3728  tran_ia_sap1_sol += tran_std_sap1_sol[k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3729  tran_ia_sap1_sen += tran_std_sap1_sen[k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3730  tran_ia_sb_sol += tran_std_sb_sol [k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3731  tran_ia_sb_sen += tran_std_sb_sen [k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3732  tran_ia_sbp1_sol += tran_std_sbp1_sol[k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3733  tran_ia_sbp1_sen += tran_std_sbp1_sen[k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
3734  // printf("%d %d TRAN_IA_SA_SOL= %g %g %g %d\n",j,k,tran_ia_sa_sol, tran_std_sa_sol [k],init_speccal15_.finstr[j][k-ia[j]+ncvhf[j]-1],ia[j]);
3735  }
3736  }
3737 
3738  ia_p1 = ia[j] + 1;
3739  for (k = ia_p1 - (ncvhf[j] - 1); k < ia_p1 + ncvhf[j]; k++) {
3740  for (i = 0; i < NH2OMAX; i++) {
3741  tran_iap1[i] += tran_std[i][k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3742  }
3743  if (splitpaths) {
3744  tran_iap1_sa_sol += tran_std_sa_sol [k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3745  tran_iap1_sa_sen += tran_std_sa_sen [k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3746  tran_iap1_sap1_sol += tran_std_sap1_sol[k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3747  tran_iap1_sap1_sen += tran_std_sap1_sen[k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3748  tran_iap1_sb_sol += tran_std_sb_sol [k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3749  tran_iap1_sb_sen += tran_std_sb_sen [k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3750  tran_iap1_sbp1_sol += tran_std_sbp1_sol[k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3751  tran_iap1_sbp1_sen += tran_std_sbp1_sen[k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
3752  // printf("%d %d TRAN_IAP1_SA_SOL= %g %g %g\n",j,k,tran_iap1_sa_sol, tran_std_sa_sol [k],init_speccal15_.finstr[j][k-ia[j]+ncvhf[j]-1]);
3753  }
3754  }
3755 
3756  // Linear interpolation to get TRNCAL from TRAN_IA and TRAN_IAP1:
3757  dlt_ia = wavln_std[ia_p1] - wavln_std[ia[j]];
3758  fia = (wavln_std[ia_p1] - wavobs[j]) / dlt_ia;
3759  fia_p1 = (float) 1.0 - fia;
3760 
3761  for (i = 0; i < NH2OMAX; i++)
3762  tran_table1_.trntbl[i][j] = fia * tran_ia[i] + fia_p1 * tran_iap1[i];
3763 
3764  if (splitpaths) {
3765  tg_sol_a = fia * tran_ia_sa_sol + fia_p1*tran_iap1_sa_sol;
3766  tg_sol_ap1 = fia * tran_ia_sap1_sol + fia_p1*tran_iap1_sap1_sol;
3767  tg_sen_a = fia * tran_ia_sa_sen + fia_p1*tran_iap1_sa_sen;
3768  tg_sen_ap1 = fia * tran_ia_sap1_sen + fia_p1*tran_iap1_sap1_sen;
3769 
3770  tg_sol_b = fia * tran_ia_sb_sol + fia_p1*tran_iap1_sb_sol;
3771  tg_sol_bp1 = fia * tran_ia_sbp1_sol + fia_p1*tran_iap1_sbp1_sol;
3772  tg_sen_b = fia * tran_ia_sb_sen + fia_p1*tran_iap1_sb_sen;
3773  tg_sen_bp1 = fia * tran_ia_sbp1_sen + fia_p1*tran_iap1_sbp1_sen;
3774 
3775  tg_sol[j] = (float) 0.5 * (f1a * tg_sol_a + f2a * tg_sol_ap1 +
3776  f1b * tg_sol_b + f2b * tg_sol_bp1);
3777  tg_sen[j] = (float) 0.5 * (f1a * tg_sen_a + f2a * tg_sen_ap1 +
3778  f1b * tg_sen_b + f2b * tg_sen_bp1);
3779  // printf("%d %d TG_SOL_SA= %g %g %g %g %g %g %g <%g %g %g %g>\n",j,ia[j],tg_sol[j],tg_sol_a,tg_sol_ap1,
3780  // tran_ia_sa_sol,tran_iap1_sa_sol,
3781  // tran_ia_sap1_sol,tran_iap1_sap1_sol,
3782  // f1a,f2a,f1b,f2b);
3783  // printf("%d %d TG_SOL_SB= %g %g %g %g %g %g %g <%g %g %g %g>\n",j,ia[j],tg_sol[j],tg_sol_b,tg_sol_bp1,
3784  // tran_ia_sb_sol,tran_iap1_sb_sol,
3785  // tran_ia_sbp1_sol,tran_iap1_sbp1_sol,
3786  // f1a,f2a,f1b,f2b);
3787  }
3788 
3789  if (!getinput5_.full_calc) {
3790  for (k = 0; k < NH2OMAX; k++) {
3791  tran_table1_.diff_tran[k][j] = (tran_table1_.trntbl[k][j] - tran_table1_.tran_kd[k][j]) /
3792  tran_table1_.tran_kd[k][j]+(float) 1.0;
3793  // printf("RJH: TRAN2: %d %d %f %f %f\n",j,k,tran_table1_.tran_kd[k][j],tran_table1_.trntbl[k][j],tran_table1_.diff_tran[k][j]);
3794  }
3795  }
3796  } //nbands
3797 
3798  for (i = 0; i < NH2OMAX; i++) {
3799  free(tran_std[i]);
3800  free(tran_med[i]);
3801  }
3802  free(tran_std);
3803  free(tran_med);
3804 }
3805 
3806 /* Converted to C from Fortran by R. Healy (richard.healy@nasa.gov) 4/2016
3807  *
3808 !********************************************************************************
3809 !* *
3810 !* Name: TRAN_SMOOTH_OTHERS *
3811 !* Purpose: This program is to smooth the line-by-line high resolution *
3812 !* spectrum to lower resolution spectrum that matches the resolutions *
3813 !* of imaging spectrometer data - for other absorbing gases *
3814 !* Parameters: none. *
3815 !* Algorithm: The smoothing is done in two stages. The 1st stage is to smooth *
3816 !* the high resolution spectrum to medium resolution spectrum at a *
3817 !* constant FWHM (0.2 nm) and a constant wavelength interval *
3818 !* (0.1 nm). The 2nd stage smoothing is to smooth the medium *
3819 !* resolution spectrum to resolutions of input imaging spectrometer *
3820 !* data. *
3821 !* Globals used: The global variables used are contained in the file *
3822 !* "COMMONS_INC" *
3823 !* Global output: TRNCAL - total transmittances of all gases that match the *
3824 !* resolutions of imaging spectrometers. *
3825 !* Return Codes: none. *
3826 !* *
3827 !********************************************************************************
3828  */
3829 void tran_smooth_others() {
3830 
3831  static int firstCall = 1;
3832 
3833  /*
3834  * Converted to C from Fortran by R. Healy (richard.healy@nasa.gov) 4/2016
3835  *
3836  First stage of smoothing - smooth line-by-line high resolution spectrum with
3837  over 300,000 points (point spacing of 0.05 cm-1) to medium resolution
3838  spectrum (resolution of 0.2 nm and point spacing of 0.1 nm) with about
3839  25,000 points.
3840 
3841  The smoothing of line-by-line spectrum is done in wavenumber domain. For
3842  a spectrum with a constant 0.2 nm resolution in wavelength domain, it has
3843  variable resolution in wavenumber domain. This effect is properly taken
3844  care of in the design of smoothing functions.
3845 
3846  Because the high resolution spectrum is in wavenumber units (cm-1), while
3847  the medium resolution spectrum is in wavelength units, the two kinds of
3848  grids do not automatically match. In order to match the grids, arrays
3849  of INDEX_MED and TRAN_MED_INDEX are specially designed. The desired
3850  medium resolution spectrum, TRAN_MED, at constant 0.1 nm point spacing
3851  and 0.2 nm resolution is obtained through linear interpolation of
3852  TRAN_MED_INDEX array.
3853 
3854  */
3855  int32_t i, j, k;
3856  float *wavln_med_index = init_speccal13_.wavln_med_index, *fwhm = getinput4_.fwhm, *wavln_std = init_speccal12_.wavln_std;
3857  float *wavobs = getinput4_.wavobs;
3858  static float **finstr_wavno, *fwhm_wavno;
3859  static int32_t *ncvhf_wavno, *ia, ia_p1;
3860  int32_t ncvtot_wavno, ncvtot, nbands = getinput5_.nbands, *ncvhf = init_speccal15_.ncvhf, *index_med = init_speccal13_.index_med;
3861 
3862  float *trntblo = tran_table1_.trntblo;
3863  float *tg_solo = tran_table_l2gen_.tg_solo, *tg_seno = tran_table_l2gen_.tg_seno;
3864 
3865  float *tran_hi_others = init_speccal1_.tran_hi_others;
3866  float tran_ia, tran_iap1;
3867  float tran_ias_sol, tran_ias_sen, tran_iap1s_sol, tran_iap1s_sen;
3868  float sumins, dlt, fj, fjm1, dlt_ia, fia, fia_p1;
3869  // Transmittance of medium resolution data
3870  float *tran_std_o, *tran_med_o, *tran_med_index_o;
3871  float *tran_std_os_sol, *tran_std_os_sen;
3872  float *tran_med_os_sol, *tran_med_os_sen, *tran_med_index_os_sol, *tran_med_index_os_sen;
3873  int splitpaths = geometry_l2gen_.splitpaths;
3874  float airmass, airm, mu = geometry5_.mu, mu0 = geometry5_.mu0;
3875  float tran_hi_others_sol, tran_hi_others_sen;
3876 
3877 
3878  int32_t np_std = NP_STD, ndx1, ndx2;
3879 
3880  if (firstCall) {
3881  if (!(fwhm_wavno = (float *) malloc(NP_MED * sizeof (float)))) {
3882  printf("-E- : Error allocating memory to fwhm_wavno\n");
3883  exit(FATAL_ERROR);
3884  }
3885  if (!(finstr_wavno = (float **) malloc(NP_MED * sizeof (float*)))) {
3886  printf("-E- : Error allocating memory to finstr_wavno\n");
3887  exit(FATAL_ERROR);
3888  }
3889  ncvhf_wavno = (int32_t*) allocateMemory(NP_MED * sizeof (int32_t), "ncvhf_wavno");
3890 
3891  for (j = 0; j < NP_MED; j++) {
3892  if (!(finstr_wavno[j] = (float *) malloc(NINSTRF * sizeof (float)))) {
3893  printf("-E- : Error allocating memory to finstr_wavno\n");
3894  exit(FATAL_ERROR);
3895  }
3896 
3897  fwhm_wavno[j] = (float) 10000. * DLT_MED / (wavln_med_index[j] * wavln_med_index[j]);
3898  ncvhf_wavno[j] = (FACDLT * fwhm_wavno[j] / DWAVNO + 1.);
3899  ncvtot_wavno = 2 * ncvhf_wavno[j] - 1;
3900 
3901  sumins = 0;
3902  for (i = ncvhf_wavno[j] - 1; i < ncvtot_wavno; i++) {
3903  finstr_wavno[j][i] = exp(-CONST1 * pow((float) (i - ncvhf_wavno[j] + 1) * DWAVNO / fwhm_wavno[j], 2));
3904  sumins += finstr_wavno[j][i];
3905  // if (j==0) printf("A:finstr_wavno sumins: %d %g %g %d\n",i,sumins,finstr_wavno[j][i],i-ncvhf_wavno[j]+1);
3906  }
3907 
3908  for (i = 0; i < ncvhf_wavno[j] - 1; i++) {
3909  finstr_wavno[j][i] = finstr_wavno[j][ncvtot_wavno - i - 1];
3910  sumins += finstr_wavno[j][i];
3911  // if (j==0) printf("B:finstr_wavno sumins: %d %g %g %d\n",i,sumins,finstr_wavno[j][i],ncvtot_wavno-i-1);
3912  }
3913 
3914  sumins *= DWAVNO;
3915 
3916  for (i = 0; i < ncvtot_wavno; i++) {
3917  finstr_wavno[j][i] *= (DWAVNO / sumins);
3918  // printf("RJH: finstr_wavno %d %d %g %g \n",i,j,finstr_wavno[j][i],sumins);
3919  }
3920 
3921  }
3922 
3923  ia = (int32_t*) allocateMemory(nbands * sizeof (int32_t), "ia");
3924  //Calculate instrumental response functions...
3925 
3926  for (j = 0; j < nbands; j++) {
3927  sumins = 0.0;
3928  ncvtot = 2 * ncvhf[j] - 1;
3929  // if (j==0) printf("NCVTOT: %d %d\n",ncvtot,ncvhf[j]);
3930  for (i = ncvhf[j] - 1; i < ncvtot; i++) {
3931  init_speccal15_.finstr[j][i] = exp(-CONST1 * pow((float) (i - ncvhf[j] + 1) * DWAVLN / fwhm[j], 2));
3932  sumins += init_speccal15_.finstr[j][i];
3933  // if (j==0) printf("c:A:finstr sumins: %d %g %g %d fwhm=%g %g %g %g %g %g\n",i,sumins,init_speccal15_.finstr[j][i],i-ncvhf[j]+1,fwhm[j],DWAVLN,
3934  // (float) (i-ncvhf[j]+1)*DWAVLN/fwhm[j],pow( (float) (i-ncvhf[j]+1)*DWAVLN/fwhm[j],2),-CONST1*pow( (float) (i-ncvhf[j]+1)*DWAVLN/fwhm[j],2),
3935  // exp(-CONST1*pow( (float) (i-ncvhf[j]+1)*DWAVLN/fwhm[j],2)));
3936  }
3937 
3938  for (i = 0; i < ncvhf[j] - 1; i++) {
3939  init_speccal15_.finstr[j][i] = init_speccal15_.finstr[j][ncvtot - i - 1];
3940  sumins += init_speccal15_.finstr[j][i];
3941  // if (j==0) printf("c:B:finstr sumins: %d %g %g %d\n",i,sumins,init_speccal15_.finstr[j][i],ncvtot-i-1);
3942  }
3943 
3944  sumins *= DWAVLN;
3945 
3946  for (i = 0; i < ncvtot; i++)
3947  init_speccal15_.finstr[j][i] *= (DWAVLN / sumins);
3948 
3949  ia[j] = hunt(wavln_std, np_std, (double) wavobs[j], ia[j]);
3950  }
3951 
3952  firstCall = 0;
3953  }
3954 
3955  airmass = 1.0 / mu0 + 1.0 / mu;
3956 
3957  /*
3958  * !!!High resolution transmittances for CO2, N2O, CO, CH4, and O2 should
3959  also be calculated for wavelength start = 0.56 micron.
3960  *
3961  */
3962 
3963  if (!(tran_std_o = (float *) calloc(NP_STD, sizeof (float)))) {
3964  printf("-E- : Error allocating memory to tran_std_o\n");
3965  exit(FATAL_ERROR);
3966  }
3967  if (!(tran_std_os_sol = (float *) calloc(NP_STD, sizeof (float)))) {
3968  printf("-E- : Error allocating memory to tran_std_os_sol\n");
3969  exit(FATAL_ERROR);
3970  }
3971  if (!(tran_std_os_sen = (float *) calloc(NP_STD, sizeof (float)))) {
3972  printf("-E- : Error allocating memory to tran_std_os_sen\n");
3973  exit(FATAL_ERROR);
3974  }
3975  if (!(tran_med_o = (float *) calloc(NP_MED, sizeof (float)))) {
3976  printf("-E- : Error allocating memory to tran_med_o\n");
3977  exit(FATAL_ERROR);
3978  }
3979  if (!(tran_med_index_o = (float *) calloc(NP_MED, sizeof (float)))) {
3980  printf("-E- : Error allocating memory to tran_med_index_o\n");
3981  exit(FATAL_ERROR);
3982  }
3983  if (!(tran_med_os_sol = (float *) calloc(NP_MED, sizeof (float)))) {
3984  printf("-E- : Error allocating memory to tran_med_os_sol\n");
3985  exit(FATAL_ERROR);
3986  }
3987  if (!(tran_med_os_sen = (float *) calloc(NP_MED, sizeof (float)))) {
3988  printf("-E- : Error allocating memory to tran_med_os_sen\n");
3989  exit(FATAL_ERROR);
3990  }
3991  if (!(tran_med_index_os_sol = (float *) calloc(NP_MED, sizeof (float)))) {
3992  printf("-E- : Error allocating memory to tran_med_index_os_sol\n");
3993  exit(FATAL_ERROR);
3994  }
3995  if (!(tran_med_index_os_sen = (float *) calloc(NP_MED, sizeof (float)))) {
3996  printf("-E- : Error allocating memory to tran_med_index_os_sen\n");
3997  exit(FATAL_ERROR);
3998  }
3999 
4000  for (j = 0; j < NP_MED; j++) {
4001  ndx1 = index_med[j]-(ncvhf_wavno[j] - 1);
4002  ndx2 = index_med[j] + ncvhf_wavno[j] - 1;
4003 
4004  //printf(" NDX1,NDX2= %d %d %d %d\n",ndx1,ndx2,index_med[j],ncvhf_wavno[j]);
4005  for (k = ndx1 - 1; k < ndx2; k++) {
4006  tran_med_index_o[j] += tran_hi_others[k] * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
4007 
4008  // printf(" TRAN_MED_INDEX_O: %d %d <%g> %g %g %d %d %d %g\n",j,k,tran_med_index_o[j],tran_hi_others[k],
4009  // finstr_wavno[j][k-index_med[j]+ncvhf_wavno[j]],
4010  // index_med[j],ncvhf_wavno[j],k-index_med[j]+ncvhf_wavno[j],tran_hi_others[k]*finstr_wavno[j][k-index_med[j]+ncvhf_wavno[j]]);
4011  if (splitpaths) {
4012  //Note that the tran_hi_* arrays are 0 for solar path and 1 for sensor
4013  airm = -(float) log(tran_hi_others[k]) / airmass;
4014  tran_hi_others_sol = (float) exp(-airm / mu0);
4015  tran_hi_others_sen = (float) exp(-airm / mu);
4016  tran_med_index_os_sol[j] += tran_hi_others_sol * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
4017  tran_med_index_os_sen[j] += tran_hi_others_sen * finstr_wavno[j][k - index_med[j] + ncvhf_wavno[j]];
4018 
4019  }
4020  }
4021  }
4022  /*
4023  *
4024  Linear interpolation to get TRAN_MED from TRAN_MED_INDEX:
4025  (Note that WAVLN_MED_INDEX(J) >= WAVLN_MED(J) )
4026  *
4027  */
4028  tran_med_o[0] = tran_med_index_o[0];
4029  tran_med_o[NP_MED - 1] = tran_med_index_o[NP_MED - 1];
4030 
4031  if (splitpaths) {
4032  tran_med_os_sol[0] = tran_med_index_os_sol[0];
4033  tran_med_os_sol[NP_MED - 1] = tran_med_index_os_sol[NP_MED - 1];
4034  tran_med_os_sen[0] = tran_med_index_os_sen[0];
4035  tran_med_os_sen[NP_MED - 1] = tran_med_index_os_sen[NP_MED - 1];
4036  }
4037 
4038  for (j = 1; j < NP_MED - 1; j++) {
4039  if (wavln_med_index[j] <= init_speccal12_.wavln_med[j]) {
4040  tran_med_o[j] = tran_med_index_o[j];
4041  // printf(" TRAN_MED_WAVLN: %d <%g> <%g> %g %g\n",j,tran_med_o[j],tran_med_index_o][j],
4042  // wavln_med_index[j],init_speccal12_.wavln_med[j]);
4043  if (splitpaths) {
4044  tran_med_os_sol[j] = tran_med_index_os_sol[j];
4045  tran_med_os_sen[j] = tran_med_index_os_sen[j];
4046  }
4047  } else {
4048  dlt = wavln_med_index[j] - wavln_med_index[j - 1];
4049  fjm1 = (wavln_med_index[j] - init_speccal12_.wavln_med[j]) / dlt;
4050  fj = (init_speccal12_.wavln_med [j] - wavln_med_index[j - 1]) / dlt;
4051  tran_med_o[j] = fjm1 * tran_med_index_o[j - 1] + fj * tran_med_index_o[j];
4052  // printf(" TRAN_MED_INTERP: %d <%g> <%g> %g %g %g\n",j,tran_med_index_o[j-1],tran_med_index_o[j],
4053  // fjm1,fj,dlt);
4054  if (splitpaths) {
4055  tran_med_os_sol[j] = fjm1 * tran_med_index_os_sol[j - 1] + fj * tran_med_index_os_sol[j];
4056  tran_med_os_sen[j] = fjm1 * tran_med_index_os_sen[j - 1] + fj * tran_med_index_os_sen[j];
4057  // printf(" TRAN_MED_INTERP_OS: %d <%g> <%g> %g %g %g\n",j,tran_med_index_os_sol[j-1],tran_med_index_os_sol[j],
4058  // fjm1,fj,dlt);
4059  }
4060  }
4061  }
4062 
4063  /*
4064  * --- Here multiplying O3 and NO2 spectra and other spectrum at medium resolution:
4065  *
4066  */
4067 
4068  for (i = 0; i < NP_STD; i++) tran_std_o[i] = 1.0;
4069 
4070  if (splitpaths) {
4071  for (i = 0; i < NP_STD; i++) {
4072  tran_std_os_sol[i] = 1.0;
4073  tran_std_os_sen[i] = 1.0;
4074  }
4075  }
4076 
4077  for (i = 0; i < NO3PT; i++) tran_std_o[i] = init_speccal16_.tran_o3_std[i] * init_speccal17_.tran_no2_std[i];
4078 
4079  if (splitpaths) {
4080  for (i = 0; i < NO3PT; i++) {
4081  airm = -(float) log(tran_std_o[i]) / airmass;
4082  tran_std_os_sol[i] = (float) exp(-airm / mu0);
4083  tran_std_os_sen[i] = (float) exp(-airm / mu);
4084  }
4085  }
4086 
4087  for (i = NPSHIF; i < NP_STD; i++) tran_std_o[i] *= tran_med_o[i - NPSHIF];
4088 
4089  if (splitpaths) {
4090  for (i = NPSHIF; i < NP_STD; i++) {
4091  tran_std_os_sol[i] *= tran_med_os_sol[i - NPSHIF];
4092  tran_std_os_sen[i] *= tran_med_os_sen[i - NPSHIF];
4093  // printf("%d TRAN_STD_OS_SOL= %g %g\n",i,tran_std_os_sol[i],tran_med_os_sol[i-NPSHIF]);
4094  // printf("%d TRAN_STD_OS_SEN= %g %g\n",i,tran_std_os_sen[i],tran_med_os_sen[i-NPSHIF]);
4095  }
4096  }
4097 
4098  /*
4099  * The 2nd stage of smoothing - smooth the medium resolution spectrum (resolution
4100  of 0.2 nm and point spacing of 0.1 nm) with about 25,000 points to match
4101  the coarser and variable resolution spectrum from imaging spectrometers.
4102 
4103  Initialize some index parameters:
4104 
4105  */
4106 
4107  for (j = 0; j < nbands; j++) {
4108 
4109  tran_ia = 0.0;
4110  tran_iap1 = 0.0;
4111 
4112  if (splitpaths) {
4113  tran_ias_sol = 0.0;
4114  tran_iap1s_sol = 0.0;
4115  tran_ias_sen = 0.0;
4116  tran_iap1s_sen = 0.0;
4117  }
4118 
4119  // Smoothing ...
4120  for (k = ia[j]-(ncvhf[j] - 1); k < ia[j] + ncvhf[j]; k++) {
4121  tran_ia += tran_std_o[k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
4122  // printf("TRAN_IA: %d %d %d %g %g %d %g %d\n",j,k,ia[j],tran_ia,tran_std_o[k],k-ia[j]+ncvhf[j]-1,init_speccal15_.finstr[j][k-ia[j]+ncvhf[j]-1],ncvhf[j]);
4123 
4124  if (splitpaths) {
4125  tran_ias_sol += tran_std_os_sol [k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
4126  tran_ias_sen += tran_std_os_sen [k] * init_speccal15_.finstr[j][k - ia[j] + ncvhf[j] - 1];
4127  // printf("TRAN_IAS_SOL: %d %d %g %g %d %g %d %d\n",j,k,tran_ias_sol,tran_std_os_sol[k],k-ia[j]+ncvhf[j]-1,mu0,ia[j],ncvhf[j]);
4128  // printf("TRAN_IAS_SEN: %d %d %g %g %d %g %d %d\n",j,k,tran_ias_sen,tran_std_os_sen[k],k-ia[j]+ncvhf[j]-1,mu,ia[j],ncvhf[j]);
4129  }
4130  }
4131 
4132  ia_p1 = ia[j] + 1;
4133  for (k = ia_p1 - (ncvhf[j] - 1); k < ia_p1 + ncvhf[j]; k++) {
4134  tran_iap1 += tran_std_o[k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
4135  if (splitpaths) {
4136  tran_iap1s_sol += tran_std_os_sol [k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
4137  tran_iap1s_sen += tran_std_os_sen [k] * init_speccal15_.finstr[j][k - ia_p1 + ncvhf[j] - 1];
4138  }
4139  }
4140 
4141  // Linear interpolation to get TRNCAL from TRAN_IA and TRAN_IAP1:
4142  dlt_ia = wavln_std[ia_p1] - wavln_std[ia[j]];
4143  fia = (wavln_std[ia_p1] - wavobs[j]) / dlt_ia;
4144  fia_p1 = (float) 1.0 - fia;
4145 
4146  trntblo[j] = fia * tran_ia + fia_p1*tran_iap1;
4147 
4148  if (splitpaths) {
4149 
4150  tg_solo[j] = fia * tran_ias_sol + fia_p1*tran_iap1s_sol;
4151  tg_seno[j] = fia * tran_ias_sen + fia_p1*tran_iap1s_sen;
4152  // printf("TG_SOLO: %d %g %g %g %g %g\n",j,tg_solo[j],tran_ias_sol,tran_iap1s_sol,fia,fia_p1);
4153  // printf("TG_SENO: %d %g %g %g %g %g\n",j,tg_seno[j],tran_ias_sen,tran_iap1s_sen,fia,fia_p1);
4154  }
4155  }
4156 
4157  free(tran_std_o);
4158  free(tran_med_o);
4159  free(tran_std_os_sol);
4160  free(tran_std_os_sen);
4161  free(tran_med_index_o);
4162  free(tran_med_os_sol);
4163  free(tran_med_os_sen);
4164  free(tran_med_index_os_sol);
4165  free(tran_med_index_os_sen);
4166 }
4167 
4168 //THE END!!!
clmvapp
float clmvapp
Definition: atrem_corl1.h:178
tran_std_sap1_sol
float tran_std_sap1_sol[NP_STD]
Definition: atrem_corl1v2.h:178
tg_seno
float tg_seno[NBANDS]
Definition: atrem_corl1.h:127
ATREM_N2O
#define ATREM_N2O
Definition: atrem_corl1.h:28
tg_sen
float tg_sen[NBANDS]
Definition: atrem_corl1.h:127
get_atrem_cor
int get_atrem_cor(l1str *l1rec, int32_t ip, float *rhot, float *tg_tot, float *tg_sol, float *tg_sen)
Definition: get_atrem_corl1v3.c:147
f1a
float f1a
Definition: atrem_corl1.h:229
init_tpvmr_nc
int init_tpvmr_nc(int model)
Definition: get_atrem_corl1v3.c:734
ecdf_
void ecdf_(float *xcdf, float *ycdf, int32_t *bin_number, float *xs, int32_t *sample_size)
Definition: numerical.c:22
f1b
float f1b
Definition: atrem_corl1.h:229
q
data_t q
Definition: decode_rs.h:74
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
sensorId2SensorDir
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:240
tran_hi_others
float tran_hi_others[NP_HI]
Definition: atrem_corl1v2.h:240
ist4
int32_t ist4
Definition: atrem_cor.h:110
t
data_t t[NROOTS+1]
Definition: decode_rs.h:77
tran_std_sbp1_sen
float tran_std_sbp1_sen[NP_STD]
Definition: atrem_corl1v2.h:180
tran_std_sb_sen
float tran_std_sb_sen[NP_STD]
Definition: atrem_corl1v2.h:179
co
int32_t co
Definition: atrem_corl1.h:118
RAD_DEG
#define RAD_DEG
Definition: atrem_corl1v2.h:48
getinput4_
struct @25 getinput4_
channelRatio
void channelRatio()
Definition: get_atrem_corl1v3.c:1650
H2O_BIT
#define H2O_BIT
Definition: l12_parms.h:60
j
int j
Definition: decode_rs.h:73
wavln_med_index
float wavln_med_index[NP_MED]
Definition: atrem_corl1v2.h:280
day
int32_t day
Definition: atrem_corl1v2.h:308
ied4
int32_t ied4
Definition: atrem_cor.h:110
dvap_layer
float dvap_layer
Definition: atrem_corl1.h:178
tran_med_index_sa_sol
float tran_med_index_sa_sol[NP_MED]
Definition: atrem_corl1v2.h:169
splitpaths
int splitpaths
Definition: atrem_corl1.h:228
CONST1
#define CONST1
Definition: atrem_corl1v2.h:49
status
int status
Definition: l1_czcs_hdf.c:32
getinput3_
struct @24 getinput3_
dp_layer
float dp_layer
Definition: atrem_corl1.h:178
input_l2gen_
struct @22 input_l2gen_
list
list(APPEND LIBS ${PGSTK_LIBRARIES}) add_executable(atteph_info_modis atteph_info_modis.c) target_link_libraries(atteph_info_modis $
Definition: CMakeLists.txt:7
go3
float go3
Definition: atrem_corl1v2.h:312
nl
int32_t nl
Definition: atrem_corl1.h:132
NO3PT
#define NO3PT
Definition: atrem_corl1v2.h:22
model_adj1_
struct @31 model_adj1_
o3
int32_t o3
Definition: atrem_corl1.h:118
ncvhf
int32_t ncvhf[NBANDS]
Definition: atrem_corl1v2.h:290
geometry1_
struct @77 geometry1_
f1
float f1(float x)
get_atrem
float get_atrem(float *tg_tot, float *rhot, paramstr *P)
Definition: get_atrem_corl1v3.c:346
get_angle_limits
int get_angle_limits(float **anglelimit, float **insenz, float **insolz, int *n_senz, int *n_solz)
Definition: get_atrem_corl1v3.c:1384
mapgen.levels
list levels
Definition: mapgen.py:180
wndow4
float wndow4
Definition: atrem_corl1.h:153
hp
float hp[MODELMAX]
Definition: atrem_corl1.h:173
co2
int32_t co2
Definition: atrem_corl1.h:118
allocateMemory
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
nh2o
int32_t nh2o
Definition: atrem_cor.h:102
totlo3
float totlo3
Definition: atrem_corl1v2.h:312
NULL
#define NULL
Definition: decode_rs.h:63
sno2
float sno2
Definition: atrem_corl1.h:133
tran_med_sbp1_sen
float tran_med_sbp1_sen[NP_MED]
Definition: atrem_corl1v2.h:176
NP_STD
#define NP_STD
Definition: atrem_corl1v2.h:27
tran_table1_
struct @9 tran_table1_
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
wndow1
float wndow1
Definition: atrem_corl1.h:153
wndow2
float wndow2
Definition: atrem_corl1.h:153
CO2_BIT
#define CO2_BIT
Definition: l12_parms.h:58
DLT_MED
#define DLT_MED
Definition: atrem_corl1v2.h:35
geometry
void geometry()
Definition: get_atrem_corl1v3.c:2147
tran_med_index_sa_sen
float tran_med_index_sa_sen[NP_MED]
Definition: atrem_corl1v2.h:169
wavln_std
float wavln_std[NP_STD]
Definition: atrem_corl1v2.h:275
wt4
float wt4
Definition: atrem_cor.h:114
init_speccal17_
struct @75 init_speccal17_
nb4
int32_t nb4
Definition: atrem_cor.h:98
tran_med_index_sb_sol
float tran_med_index_sb_sol[NP_MED]
Definition: atrem_corl1v2.h:171
tran_smooth
void tran_smooth(float *tran_hi)
Definition: get_atrem_corl1v3.c:3375
mu
float mu
Definition: atrem_corl1v2.h:324
init_speccal5_
struct @65 init_speccal5_
h
float h[MODELMAX]
Definition: atrem_corl1.h:131
get_tpvmr
void get_tpvmr(size_t layers, size_t models, int sds_id, char filename[FILENAME_MAX], char *varname, float *var_a)
Definition: get_atrem_corl1v3.c:707
lat
float * lat
Definition: l1_hmodis_hdf.c:756
iedp94
int32_t iedp94
Definition: atrem_cor.h:106
tm
float tm[MODELMAX]
Definition: atrem_corl1v2.h:248
O3_BIT
#define O3_BIT
Definition: l12_parms.h:57
vaptot
float vaptot[TBLMAX]
Definition: atrem_cor.h:128
nbp94
int32_t nbp94
Definition: atrem_cor.h:98
tran_med_index_sap1_sen
float tran_med_index_sap1_sen[NP_MED]
Definition: atrem_corl1v2.h:170
vmrm
float vmrm[MODELMAX]
Definition: atrem_corl1v2.h:248
debug_atrem
struct @46 debug_atrem
tran_med_sa_sol
float tran_med_sa_sol[NP_MED]
Definition: atrem_corl1v2.h:173
tran_med_sb_sol
float tran_med_sb_sol[NP_MED]
Definition: atrem_corl1v2.h:175
fwhm
float fwhm[NBANDS]
Definition: atrem_corl1.h:144
trntblo
float trntblo[NBANDS]
Definition: atrem_corl1v2.h:303
getinput14_
struct @30 getinput14_
init_speccal6_
struct @4 init_speccal6_
vrto3
float vrto3
Definition: atrem_corl1.h:133
clmvap
float clmvap
Definition: atrem_corl1.h:169
init_spectral_calculations
void init_spectral_calculations()
Definition: get_atrem_corl1v3.c:2436
hunt
int32_t hunt(float *xx, int32_t n, double x, int32_t jlo)
Definition: get_atrem_corl1v3.c:619
tp
float tp[MODELMAX]
Definition: atrem_corl1.h:173
get_abscf_data
void get_abscf_data(int levels, int bands, int sds_id, char filename[FILENAME_MAX], float *abscf, char *varname)
Definition: get_atrem_corl1v3.c:2296
findMatch
int32_t findMatch(float *list, int32_t nbands, float elem)
Definition: get_atrem_corl1v3.c:1605
tran_med_sbp1_sol
float tran_med_sbp1_sol[NP_MED]
Definition: atrem_corl1v2.h:176
spacetrack_tles.year
year
Definition: spacetrack_tles.py:170
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
tran_med_index_sbp1_sen
float tran_med_index_sbp1_sen[NP_MED]
Definition: atrem_corl1v2.h:172
hico.proc_hico.senz
senz
Definition: proc_hico.py:239
input
instr * input
Definition: main_l2binmatch.cpp:27
names::if
character(len=1000) if
Definition: names.f90:13
MODELMAX
#define MODELMAX
Definition: atrem_corl1.h:20
h2o
int32_t h2o
Definition: atrem_corl1.h:118
rno2cf
float rno2cf[NO2PT]
Definition: atrem_corl1v2.h:344
geometry5_
struct @81 geometry5_
jb
int jb
Definition: atrem_corl1.h:227
f2
float f2(float y)
dp_plane
float dp_plane
Definition: atrem_corl1.h:178
atrem_corl1v3.h
init_tpvmr
int init_tpvmr(int model)
Definition: get_atrem_corl1v3.c:672
tran_std_sa_sen
float tran_std_sa_sen[NP_STD]
Definition: atrem_corl1v2.h:177
bindex_get
int bindex_get(int32_t wave)
Definition: windex.c:45
dlt2
float dlt2
Definition: atrem_cor.h:94
ggeom
float ggeom
Definition: atrem_corl1v2.h:312
NPSHIF
#define NPSHIF
Definition: atrem_corl1v2.h:29
init_speccal8_
struct @6 init_speccal8_
NINSTRF
#define NINSTRF
Definition: atrem_corl1v2.h:24
no2
int32_t no2
Definition: atrem_corl1.h:118
ATREM_NO2
#define ATREM_NO2
Definition: atrem_corl1.h:24
f2a
float f2a
Definition: atrem_corl1.h:229
linterp
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
NH2OMAX
#define NH2OMAX
Definition: atrem_corl1.h:17
getinput1_
struct @21 getinput1_
nb1
int32_t nb1
Definition: atrem_cor.h:98
model_adj3_
struct @33 model_adj3_
check_polarwander.infile
infile
Definition: check_polarwander.py:71
hico.value_locate.x
list x
Definition: value_locate.py:64
init_speccal10_
struct @7 init_speccal10_
tran_table
void tran_table()
Definition: get_atrem_corl1v3.c:3019
tran_std_sbp1_sol
float tran_std_sbp1_sol[NP_STD]
Definition: atrem_corl1v2.h:180
ATREM_O3
#define ATREM_O3
Definition: atrem_corl1.h:22
nb2
int32_t nb2
Definition: atrem_cor.h:98
ATREM_CO
#define ATREM_CO
Definition: atrem_corl1.h:25
vmrp
float vmrp[MODELMAX]
Definition: atrem_corl1.h:173
getinput5_
struct @1 getinput5_
ied1p14
int32_t ied1p14
Definition: atrem_cor.h:110
FACDLT
#define FACDLT
Definition: atrem_corl1v2.h:36
tg_sol
float tg_sol[NBANDS]
Definition: atrem_corl1.h:127
pp
float pp[MODELMAX]
Definition: atrem_corl1.h:173
tran_std_sa_sol
float tran_std_sa_sol[NP_STD]
Definition: atrem_corl1v2.h:177
model_adj4_
struct @34 model_adj4_
tran_tables_
struct @50 tran_tables_
color_dtdb.nbins
int nbins
Definition: color_dtdb.py:277
getModelNum
int getModelNum(float lat, int32_t day)
Definition: get_atrem_corl1v3.c:837
NO2_BIT
#define NO2_BIT
Definition: l12_parms.h:59
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
w1p14c
float w1p14c
Definition: atrem_corl1.h:153
solzni
float solzni
Definition: atrem_corl1v2.h:307
delta
const double delta
Definition: VcstCmnConsts.h:523
istp94
int32_t istp94
Definition: atrem_cor.h:106
DWAVNO
#define DWAVNO
Definition: atrem_corl1v2.h:34
want_verbose
int want_verbose
Definition: genutils_globals.c:3
dp
float dp[MODELMAX]
Definition: atrem_corl1v2.h:248
tran_med_sap1_sen
float tran_med_sap1_sen[NP_MED]
Definition: atrem_corl1v2.h:174
unix2yds
void unix2yds(double usec, short *year, short *day, double *secs)
filename
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
tran_med_index_sap1_sol
float tran_med_index_sap1_sol[NP_MED]
Definition: atrem_corl1v2.h:170
o2
int32_t o2
Definition: atrem_corl1.h:118
ja
float ja
Definition: atrem_cor.h:114
ssh2o
float ssh2o[NH2OMAX]
Definition: atrem_corl1v2.h:312
tran_med_sap1_sol
float tran_med_sap1_sol[NP_MED]
Definition: atrem_corl1v2.h:174
ATREM_CH4
#define ATREM_CH4
Definition: atrem_corl1.h:26
NP_HI
#define NP_HI
Definition: atrem_corl1v2.h:25
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
geometry4_
struct @43 geometry4_
ist1
int32_t ist1
Definition: atrem_cor.h:106
wt2
float wt2
Definition: atrem_cor.h:114
wp94c
float wp94c
Definition: atrem_corl1.h:153
tran_std_sb_sol
float tran_std_sb_sol[NP_STD]
Definition: atrem_corl1v2.h:179
dvap_plane
float dvap_plane
Definition: atrem_corl1.h:178
init_speccal16_
struct @74 init_speccal16_
index_med
int32_t index_med[NP_MED]
Definition: atrem_corl1v2.h:279
pm
float pm[MODELMAX]
Definition: atrem_corl1v2.h:248
apply_gas_trans
void apply_gas_trans(int32_t k_surf, int levels, int bands, float *abscf, float *dp, float *vmrm, float *tran_hi)
Definition: get_atrem_corl1v3.c:2325
finst2
float finst2[FINSTMAX]
Definition: atrem_cor.h:120
tran_tables1_
struct @51 tran_tables1_
ATREM_O2
#define ATREM_O2
Definition: atrem_corl1.h:27
geometry2_
struct @78 geometry2_
geometry_l2gen_
struct @44 geometry_l2gen_
init_speccal1_
struct @63 init_speccal1_
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
getinput6_
struct @27 getinput6_
tran_smooth_others
void tran_smooth_others()
Definition: get_atrem_corl1v3.c:3829
k_plane
int32_t k_plane
Definition: atrem_corl1.h:177
fabs
#define fabs(a)
Definition: misc.h:93
vmr
float vmr[MODELMAX]
Definition: atrem_corl1.h:131
ist2
int32_t ist2
Definition: atrem_cor.h:106
obszni
float obszni
Definition: atrem_corl1v2.h:307
init_atrem
int init_atrem(int32_t sensorID, paramstr *P, l1str *l1rec, int32_t nbands)
Definition: get_atrem_corl1v3.c:877
xppp
float xppp
Definition: atrem_corl1.h:165
wt3
float wt3
Definition: atrem_cor.h:114
mu0
float mu0
Definition: atrem_corl1v2.h:324
hsurf
float hsurf
Definition: atrem_cor.h:94
ied1
int32_t ied1
Definition: atrem_cor.h:106
geometry3_
struct @10 geometry3_
tran_std_sap1_sen
float tran_std_sap1_sen[NP_STD]
Definition: atrem_corl1v2.h:178
ch4
int32_t ch4
Definition: atrem_corl1.h:118
init_speccal12_
struct @71 init_speccal12_
wavobs
float wavobs[NBANDS]
Definition: atrem_corl1.h:144
tpvmr_init1_
struct @45 tpvmr_init1_
kdistgasabs
void kdistgasabs(float *kcdf, float *abscf, float *waveno, float *wavobs, int32_t nphi, int32_t n_layers, int32_t nbands)
Definition: get_atrem_corl1v3.c:1706
locate_pos_
void locate_pos_(float *xx, int32_t *n1, float *x1, int32_t *jj)
Definition: get_atrem_corl1v3.c:2091
tran_table_l2gen_
struct @23 tran_table_l2gen_
wndow3
float wndow3
Definition: atrem_corl1.h:153
rdatreminfo
int32_t rdatreminfo(int32_t sensorID, int32_t evalmask, const char *pname, void **pval)
Definition: rdatreminfo.c:38
create_images.bands
bands
Definition: create_images.py:91
tran_med_sa_sen
float tran_med_sa_sen[NP_MED]
Definition: atrem_corl1v2.h:173
model
int32_t model
Definition: atrem_corl1.h:132
VSTART
#define VSTART
Definition: atrem_corl1v2.h:31
for
for(i=0;i< NROOTS;i++) s[i]
Definition: decode_rs.h:85
get_current_angle_limit
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)
Definition: get_atrem_corl1v3.c:1553
dlt
float dlt
Definition: atrem_cor.h:94
NINSTR_MAX
#define NINSTR_MAX
Definition: atrem_corl1v2.h:28
getinput8_
struct @29 getinput8_
init_speccal3_
struct @3 init_speccal3_
ied3
int32_t ied3
Definition: atrem_cor.h:110
tran_med_index_sb_sen
float tran_med_index_sb_sen[NP_MED]
Definition: atrem_corl1v2.h:171
ATREM_CO2
#define ATREM_CO2
Definition: atrem_corl1.h:23
g_vap
float g_vap
Definition: atrem_cor.h:132
DWAVLN
#define DWAVLN
Definition: atrem_corl1v2.h:33
i
int i
Definition: decode_rs.h:71
model_adjust
void model_adjust()
Definition: get_atrem_corl1v3.c:1931
init_speccal11_
struct @8 init_speccal11_
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
wt1
float wt1
Definition: atrem_cor.h:114
sensorID
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:97
nb1p14
int32_t nb1p14
Definition: atrem_cor.h:98
tran_med_sb_sen
float tran_med_sb_sen[NP_MED]
Definition: atrem_corl1v2.h:175
ied2
int32_t ied2
Definition: atrem_cor.h:106
init_speccal7_
struct @5 init_speccal7_
k_surf
int32_t k_surf
Definition: atrem_corl1.h:182
NO2PT
#define NO2PT
Definition: atrem_corl1v2.h:23
init_speccal13_
struct @72 init_speccal13_
ist1p14
int32_t ist1p14
Definition: atrem_cor.h:110
k
int k
Definition: decode_rs.h:73
init_speccal15_
struct @73 init_speccal15_
NH2OMAXM1
#define NH2OMAXM1
Definition: atrem_corl1.h:18
nb3
int32_t nb3
Definition: atrem_cor.h:98
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
tran_med_index_sbp1_sol
float tran_med_index_sbp1_sol[NP_MED]
Definition: atrem_corl1v2.h:172
tg_solo
float tg_solo[NBANDS]
Definition: atrem_corl1.h:127
getinput7_
struct @2 getinput7_
o3cf
float o3cf[NO3PT]
Definition: atrem_corl1v2.h:340
ist3
int32_t ist3
Definition: atrem_cor.h:110
f2b
float f2b
Definition: atrem_corl1.h:229
nb
int32_t nb
Definition: atrem_corl1.h:132
n2o
int32_t n2o
Definition: atrem_corl1.h:118
obsphi
float obsphi
Definition: atrem_corl1v2.h:307
NP_MED
#define NP_MED
Definition: atrem_corl1v2.h:26
count
int count
Definition: decode_rs.h:79


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster