OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
par_utils.c
Go to the documentation of this file.
1 /******************************************************************************
2  *
3  * Photosynthetic Available Radiation Utilities
4  * Sensor-Independent Shared Subroutine
5  * by Robert J. Lossing
6  * July 30, 2013
7  * Ocean Color Processing Group
8  *
9  *******************************************************************************/
10 
11 /******************************************************************************
12  *
13  * Retrieve Phase function and Omega, given input scattering angle
14  * and Angstrom coefficient.
15  *
16  * Authors: Robert Frouin <RFrouin@ucsd.edu>, scientific algorithms
17  * John McPherson <JMcPherson@ucsd.edu>, program structure
18  *
19  * Note: Angstrom exponents in tables aren't in increasing order,
20  * need to determine proper order to use:
21  * Model omega(6) alpha(6,8) true_order
22  * 1 1.0000000E+00 -8.9877717E-02 2nd
23  * 2 9.9999797E-01 -9.1842167E-02 1st *min
24  * 3 9.8501903E-01 5.1216149E-01 8
25  * 4 9.8840600E-01 4.0654629E-01 7
26  * 5 9.9607497E-01 1.9749035E-01 4
27  * 6 9.9880499E-01 8.1851527E-02 3rd
28  * 7 9.7776997E-01 7.7916741E-01 9
29  * 8 9.9351001E-01 3.9742991E-01 6
30  * 9 9.9790698E-01 2.1781628E-01 5
31  * 10 9.5734900E-01 1.6471386E+00 12th *max
32  * 11 9.8160601E-01 1.4880587E+00 11th
33  * 12 9.9180502E-01 1.2930322E+00 10th
34  *
35  *******************************************************************************/
36 
37 #include <stdio.h>
38 #include <math.h>
39 #include <stdlib.h>
40 #include "par_utils.h"
41 #define NMODEL 12
42 #define NS 15
43 #define NT 9
44 
45 static float *tablewavelengths;
46 static float *tablephaseangles;
47 static float *tablealphas; //angstrom exponent
48 static float *tableomegas; //single scattering albedo
49 static float *tableaerphasefunc;
50 
51 void GetAerPhase(l2str *l2rec, int ip, int32_t nbands, float angstrom,
52  float *aerphasefunc, float *omega, float *modelAngstrom) {
53 
54  float phaseangle = l2rec->l1rec->scattang[ip];
55  float distx, dista, p1, p2;
56  int i, j, ix1=0, ix2, ia1=0, ia2, widx;
57 
58  static int firstCall = TRUE;
59 
60  // The data table is not in monotonically-increasing order, so here are the indices in proper order:
61 
62  int ii[12] = {1, 0, 5, 4, 8, 7, 3, 2, 6, 11, 10, 9};
63 
64  // Read the aerosol data file:
65 
66  if (firstCall) {
67  firstCall = FALSE;
68  tablewavelengths = (float*) allocateMemory(nbands * sizeof (float),
69  "tablewavelengths");
70  tablephaseangles = (float*) allocateMemory(NPHASE * sizeof (float),
71  "tablephaseangles");
72  tablealphas = (float*) allocateMemory(nbands * NMODEL * sizeof (float),
73  "tablealphas");
74  tableomegas = (float*) allocateMemory(nbands * NMODEL * sizeof (float),
75  "tableomegas");
76  tableaerphasefunc = (float*) allocateMemory(
77  NPHASE * nbands * NMODEL * sizeof (float), "tableaerphasefunc");
78  read_aerosol_par(l2rec, nbands, tablewavelengths, tablephaseangles,
79  tablealphas, tableomegas, tableaerphasefunc);
80  }
81 
82  // Find out which tablephaseangles-values in the table will be used:
83 
84  if (phaseangle < tablephaseangles[0]) {
85  // printf("Input phase angle too small, setting to 0. Was %f", phaseangle);
86  phaseangle = tablephaseangles[0];
87  ix1 = 0;
88  } else if (phaseangle >= tablephaseangles[NPHASE - 1]) {
89  // printf("Input phase angle too large, setting to 180. Was %f",
90  // phaseangle);
91  phaseangle = tablephaseangles[NPHASE - 1];
92  ix1 = NPHASE - 2;
93  } else {
94  for (i = NPHASE - 2; i > 0; i--) {
95  if (phaseangle >= tablephaseangles[i]) {
96  ix1 = i;
97  break;
98  }
99  }
100  }
101  ix2 = ix1 + 1;
102 
103  distx = (tablephaseangles[ix2] - phaseangle)
104  / (tablephaseangles[ix2] - tablephaseangles[ix1]);
105 
106  // Find out which models in the table to use:
107 
108  widx = windex(670.0, tablewavelengths, nbands);
109 
110  if (angstrom < tablealphas[widx * NMODEL + ii[0]]) {
111  angstrom = tablealphas[widx * NMODEL + ii[0]];
112  ia1 = 0;
113  } else if (angstrom >= tablealphas[widx * NMODEL + ii[NMODEL - 1]]) {
114  angstrom = tablealphas[widx * NMODEL + ii[NMODEL - 1]];
115  ia1 = NMODEL - 2;
116  } else {
117  for (j = NMODEL - 2; j > 0; j--) {
118  if (angstrom >= tablealphas[widx * NMODEL + ii[j]]) {
119  ia1 = j;
120  break;
121  }
122  }
123  }
124  ia2 = ia1 + 1;
125  dista = (tablealphas[widx * NMODEL + ii[ia2]] - angstrom)
126  / (tablealphas[widx * NMODEL + ii[ia2]]
127  - tablealphas[widx * NMODEL + ii[ia1]]);
128 
129  // Loop through the wavelengths and perform interpolations:
130 
131  for (i = 0; i < nbands; i++) {
132  omega[i] = (tableomegas[i * NMODEL + ii[ia1]] * dista)
133  + (tableomegas[i * NMODEL + ii[ia2]] * (1.0 - dista));
134 
135  modelAngstrom[i] = (tablealphas[i * NMODEL + ii[ia1]] * dista)
136  + (tablealphas[i * NMODEL + ii[ia2]] * (1.0 - dista));
137 
138  p1 =
139  (tableaerphasefunc[i * NMODEL * NPHASE + ii[ia1] * NPHASE + ix1]
140  * distx)
141  + (tableaerphasefunc[i * NMODEL * NPHASE
142  + ii[ia1] * NPHASE + ix2] * (1.0 - distx));
143 
144  p2 =
145  (tableaerphasefunc[i * NMODEL * NPHASE + ii[ia2] * NPHASE + ix1]
146  * distx)
147  + (tableaerphasefunc[i * NMODEL * NPHASE
148  + ii[ia2] * NPHASE + ix2] * (1.0 - distx));
149 
150  aerphasefunc[i] = (p1 * dista) + (p2 * (1.0 - dista));
151  }
152 
153  return;
154 }
155 
156 /********************************************************************************
157  *
158  * Description:
159  * This is to read the SeaWiFS aerosol data file for the PAR.
160  * Menghua Wang 6/22/99
161  * Parameters:
162  * angl, R, ---- scattering angles.
163  * omega0, R, --- single scattering albedo.
164  * alpha, R, --- Angstrom coefficient alpha(lambda,865).
165  * s11, R, --- aerosol phase function.
166  *
167  * SeaWiFS aerosol models:
168  * 1-2: Oceanic RH=90 and 99%
169  * 3-6: Maritime RH=50%, 70, 90, and 99%
170  * 7-9: Coastal RH=50, 90, and 99%
171  * 10-12: Tropospheric RH=50, 90, and 99%
172  *
173  *******************************************************************************/
174 void read_aerosol_par(l2str *l2rec, int32_t nbands, float *tablewavelengths,
175  float *tablephaseangles, float *tablealphas, float *tableomegas,
176  float *tableaerphasefunc) {
177 
178  int imodel[NMODEL];
179 
180  char aerosolfilename[FILENAME_MAX] = "";
181  char *dataroot;
182  FILE *aerosoldata;
183  //Get OCDATAROOT path from user:
184  if ((dataroot = getenv("OCDATAROOT")) == NULL) {
185  printf("OCDATAROOT environment variable is not defined.\n");
186  return;
187  }
188  sprintf(aerosolfilename, "%s/%s/%s_aerosol_par.dat", dataroot,
189  sensorId2SensorDir(l2rec->l1rec->l1file->sensorID),
190  sensorId2SensorDir(l2rec->l1rec->l1file->sensorID));
191 
192  // Opens the aerosol data file for reading.
193  if ((aerosoldata = fopen(aerosolfilename, "r")) == NULL) {
194  printf("Error reading aerosol table for PAR from %s.\n", aerosolfilename);
195  exit(EXIT_FAILURE);
196  }
197  printf("Loading aerosol properties for PAR from %s.\n", aerosolfilename);
198 
199  int j, jj, iph;
200  int num_model;
201  for (jj = 0; jj < NPHASE; jj++) {
202  if ((fscanf(aerosoldata, "%f", &tablephaseangles[jj])) == 0) {
203  printf("Problem reading phase angles from %s.\n", aerosolfilename);
204  exit(EXIT_FAILURE);
205  }
206  }
207  for (j = 0; j < nbands; j++) {
208  for (iph = 0; iph < NMODEL; iph++) {
209  if ((fscanf(aerosoldata, "%d %d %f", &imodel[iph], &num_model,
210  &tablewavelengths[j])) == 0) {
211  printf("Problem reading model number and wavelength from %s.\n",
212  aerosolfilename);
213  exit(EXIT_FAILURE);
214  }
215  if ((fscanf(aerosoldata, "%f %f", &tableomegas[j * NMODEL + iph],
216  &tablealphas[j * NMODEL + iph])) == 0) {
217  printf("Problem reading SSA and Angstrom from %s.\n",
218  aerosolfilename);
219  exit(EXIT_FAILURE);
220  }
221  for (jj = 0; jj < NPHASE; jj++) {
222  if ((fscanf(aerosoldata, "%f",
223  &tableaerphasefunc[j * NMODEL * NPHASE + iph * NPHASE
224  + jj])) == 0) {
225  printf("Problem reading phase function values from %s.\n",
226  aerosolfilename);
227  exit(EXIT_FAILURE);
228  }
229  }
230  }
231 
232  }
233  fclose(aerosoldata);
234 }
235 
236 /********************************************************************************
237  *
238  * Estimate Dobson units from climatology, given the month and
239  * latitude. Table has 12 columns, one per month, and 35 rows,
240  * from 85 N to -85, in steps of 5 deg lat.
241  * Perform 2d interpolation (latitude,month_day) to get estimate.
242  *
243  * Author: John McPherson <JMcPherson@ucsd.edu>
244  *
245  * Dobson look-up table from Keating et al., Ozone reference models
246  * for the middle atmosphere, Handbook for MAP, Vol. 31, 1-36, 1989
247  *
248  */
249 
250 float EstimateDobson(int32_t year, int32_t month, int32_t day, float lat) {
251 
252  int i1, i2, m1, m2;
253  float daymid, d1, d2;
254  float t, t1, t2, fac, fac2, difflat;
255  float dobson;
256 
257  int tabdobson[12][35] = {
258  { 395, 395, 395, 395, 395, 392, 390, 387, 376,
259  354, 322, 292, 269, 254, 248, 246, 247, 251, 255, 260, 266, 271,
260  277, 286, 295, 306, 319, 334, 344, 344, 338, 331, 324, 320, 316},
261 
262  { 433, 433, 433, 436, 432, 428, 426, 418, 402, 374, 338, 303, 278, 261, 251,
263  246, 248, 250, 254, 258, 262, 265, 270, 278, 286, 294, 303, 313,
264  322, 325, 324, 317, 306, 299, 294},
265 
266  { 467, 470, 460, 459, 451, 441, 433, 420, 401, 377, 347, 316, 291, 271, 260,
267  254, 254, 255, 257, 259, 261, 264, 269, 277, 284, 289, 296, 305,
268  312, 315, 317, 312, 305, 299, 295},
269 
270  { 467, 465, 462, 455, 444, 431, 421, 410, 395, 373, 348, 325, 304, 287, 275,
271  267, 261, 259, 258, 259, 260, 263, 271, 278, 284, 289, 297, 306,
272  314, 318, 319, 313, 302, 302, 302},
273 
274  { 411, 414, 416, 415, 410, 406, 402, 394, 382, 363, 342, 324, 307, 291, 279,
275  271, 264, 260, 258, 257, 258, 264, 271, 281, 291, 303, 312, 318,
276  322, 323, 322, 322, 322, 322, 322},
277 
278  { 371, 371, 370, 368, 367, 372, 375, 372, 360, 341, 323, 311, 301, 290, 282,
279  275, 268, 263, 259, 256, 258, 264, 273, 289, 306, 319, 327, 328,
280  328, 337, 337, 337, 337, 337, 337},
281 
282  { 333, 332, 332, 334, 338, 346, 350, 346, 335, 321, 310, 302, 296, 289, 284,
283  280, 274, 268, 262, 259, 261, 268, 279, 295, 315, 331, 340, 342,
284  338, 344, 340, 340, 340, 340, 340},
285 
286  { 311, 308, 308, 313, 320, 327, 330, 326, 319, 310, 303, 298, 291, 286, 283,
287  281, 277, 273, 268, 264, 266, 274, 288, 306, 327, 343, 353, 355,
288  351, 339, 325, 307, 294, 294, 294},
289 
290  { 283, 291, 302, 308, 312, 317, 318, 313, 307, 300, 295, 290, 284, 279, 279,
291  279, 278, 276, 272, 270, 273, 282, 295, 313, 333, 348, 360, 367,
292  368, 353, 324, 291, 267, 253, 230},
293 
294  { 299, 299, 299, 309, 315, 317, 317, 312, 302, 291, 283, 280, 275, 270, 268,
295  267, 263, 263, 265, 269, 277, 287, 301, 317, 336, 354, 371, 387,
296  402, 402, 374, 333, 294, 274, 259},
297 
298  { 314, 314, 314, 314, 332, 332, 327, 322, 311, 297, 284, 276, 270, 263, 261,
299  260, 258, 259, 264, 270, 278, 286, 298, 311, 323, 335, 350, 366,
300  381, 390, 388, 376, 357, 346, 341},
301 
302  { 358, 358, 358, 358, 358, 358, 353, 349, 338, 320, 299, 281, 267, 256, 252,
303  251, 251, 253, 257, 264, 272, 279, 287, 297, 307, 318, 332, 347,
304  358, 365, 366, 364, 358, 356, 353}
305  };
306 
307  int days[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
308 
309  // Set up for the time interpolation:
310 
311  if (year % 4 == 0)
312  days[1] = days[1] + 1;
313 
314  daymid = days[month - 1] / 2.0;
315 
316  if (day >= daymid) {
317  m1 = month - 1;
318  m2 = m1 + 1;
319  if (m2 > 12)
320  m2 = 1;
321  t = day;
322  } else {
323  m2 = month - 1;
324  m1 = m2 - 1;
325  if (m1 < 1)
326  m1 = 12;
327  t = days[m1] + day;
328  }
329 
330  t1 = days[m1] / 2.0;
331  t2 = days[m1] + (days[m2] / 2.0);
332 
333  // Perform the lat. interpolation, and further set up for time interp.:
334 
335  if (lat >= 85.0) {
336  d1 = tabdobson[m1][1];
337  d2 = tabdobson[m2][1];
338  } else if (lat <= -85.0) {
339  d1 = tabdobson[m1][34];
340  d2 = tabdobson[m2][34];
341  } else {
342  i1 = 17 - (int) (lat / 5.0);
343  i2 = i1 + 1;
344  fac = (tabdobson[month][i2] - tabdobson[month][i1]) / (-5.0);
345  difflat = lat - (90.0 - (i1 * 5.0));
346  d1 = tabdobson[m1][i1] + (fac * difflat);
347  d2 = tabdobson[m2][i1] + (fac * difflat);
348  }
349 
350  // Complete the time interpolation to get final estimate:
351 
352  fac2 = (d2 - d1) / (t2 - t1);
353 
354  dobson = d1 + (t - t1) * fac2;
355  return dobson;
356 }
357 
358 /* ****************************************************************************
359  *
360  * Interpolate water vapor table, as for dobson (ozone).
361  * Table has 12 columns, one per month, and 35 rows,
362  * from 85 N to 85 S, in steps of 5 deg lat.
363  * winter peak (Jan 15), summer peak (July 15), north hemis
364  * Perform 2d interpolation (latitude,month_day) to get estimate.
365  *
366  * Tropical: if ( ABS( Lat ) .LE. 20.0 ) WatVap = 4.12
367  *
368  * Mid-latitude case: if ( ABS( Lat ) .LE. 60.0 ) Then
369  * if ( (North.And.Summer) .Or. (South.And.Winter) ) WatVap = 2.93
370  * else WatVap = 0.85
371  *
372  * Sub-arctic case:
373  * if ( (North.And.Summer) .Or. (South.And.Winter) ) WatVap = 2.10
374  * else WatVap = 0.42
375  *
376  * after McClatchey, R.A., et al, _Optical Properties of the Atmosphere_,
377  * AFCRL 71-0279 354, 91 pp., 1971
378  *
379  */
380 
381 float EstimateWatVap(int32_t year, int32_t month, int32_t day, float lat) {
382 
383  int i1, i2, m1, m2;
384  month -= 1;
385  float daymid, d1, d2;
386  float t, t1, t2, fac, fac2, difflat;
387  float watvap;
388 
389  float tabwv[12][35] = {
390  { 0.42, 0.47, 0.52, 0.56, 0.61, 0.66, 0.71, 0.75,
391  0.80, 0.85, 1.26, 1.67, 2.08, 2.48, 2.89, 3.30, 3.71, 4.12, 3.97,
392  3.82, 3.67, 3.53, 3.38, 3.23, 3.08, 2.93, 2.84, 2.75, 2.65, 2.56,
393  2.47, 2.38, 2.28, 2.19, 2.10},
394 
395  { 0.70, 0.76, 0.81, 0.87, 0.92, 0.98, 1.03, 1.09, 1.14, 1.20, 1.56, 1.93,
396  2.29, 2.66, 3.02, 3.39, 3.75, 4.12, 3.93, 3.74, 3.54, 3.35, 3.16,
397  2.97, 2.78, 2.58, 2.50, 2.41, 2.33, 2.24, 2.16, 2.07, 1.99, 1.90,
398  1.82},
399 
400  { 0.98, 1.04, 1.11, 1.17, 1.23, 1.29, 1.36, 1.42, 1.48, 1.54, 1.87, 2.19,
401  2.51, 2.83, 3.15, 3.48, 3.80, 4.12, 3.88, 3.65, 3.41, 3.18, 2.94,
402  2.71, 2.47, 2.24, 2.16, 2.08, 2.00, 1.93, 1.85, 1.77, 1.69, 1.62,
403  1.54},
404 
405  { 1.26, 1.33, 1.40, 1.47, 1.54, 1.61, 1.68, 1.75, 1.82, 1.89, 2.17, 2.45,
406  2.73, 3.01, 3.28, 3.56, 3.84, 4.12, 3.84, 3.56, 3.28, 3.01, 2.73,
407  2.45, 2.17, 1.89, 1.82, 1.75, 1.68, 1.61, 1.54, 1.47, 1.40, 1.33,
408  1.26},
409 
410  { 1.54, 1.62, 1.69, 1.77, 1.85, 1.93, 2.00, 2.08, 2.16, 2.24, 2.47, 2.71,
411  2.94, 3.18, 3.41, 3.65, 3.88, 4.12, 3.80, 3.48, 3.15, 2.83, 2.51,
412  2.19, 1.87, 1.54, 1.48, 1.42, 1.36, 1.29, 1.23, 1.17, 1.11, 1.04,
413  0.98},
414 
415  { 1.82, 1.90, 1.99, 2.07, 2.16, 2.24, 2.33, 2.41, 2.50, 2.58, 2.78, 2.97,
416  3.16, 3.35, 3.54, 3.74, 3.93, 4.12, 3.75, 3.39, 3.02, 2.66, 2.29,
417  1.93, 1.56, 1.20, 1.14, 1.09, 1.03, 0.98, 0.92, 0.87, 0.81, 0.76,
418  0.70},
419 
420  { 2.10, 2.19, 2.28, 2.38, 2.47, 2.56, 2.65, 2.75, 2.84, 2.93, 3.08, 3.23,
421  3.38, 3.53, 3.67, 3.82, 3.97, 4.12, 3.71, 3.30, 2.89, 2.49, 2.08,
422  1.67, 1.26, 0.85, 0.80, 0.75, 0.71, 0.66, 0.61, 0.56, 0.52, 0.47,
423  0.42},
424 
425  { 1.82, 1.90, 1.99, 2.07, 2.16, 2.24, 2.33, 2.41, 2.50, 2.58, 2.78, 2.97,
426  3.16, 3.35, 3.54, 3.74, 3.93, 4.12, 3.75, 3.39, 3.02, 2.66, 2.29,
427  1.93, 1.56, 1.20, 1.14, 1.09, 1.03, 0.98, 0.92, 0.87, 0.81, 0.76,
428  0.70},
429 
430  { 1.54, 1.62, 1.69, 1.77, 1.85, 1.93, 2.00, 2.08, 2.16, 2.24, 2.47, 2.71,
431  2.94, 3.18, 3.41, 3.65, 3.88, 4.12, 3.80, 3.48, 3.15, 2.83, 2.51,
432  2.19, 1.87, 1.54, 1.48, 1.42, 1.36, 1.29, 1.23, 1.17, 1.11, 1.04,
433  0.98},
434 
435  { 1.26, 1.33, 1.40, 1.47, 1.54, 1.61, 1.68, 1.75, 1.82, 1.89, 2.17, 2.45,
436  2.73, 3.01, 3.28, 3.56, 3.84, 4.12, 3.84, 3.56, 3.28, 3.01, 2.73,
437  2.45, 2.17, 1.89, 1.82, 1.75, 1.68, 1.61, 1.54, 1.47, 1.40, 1.33,
438  1.26},
439 
440  { 0.98, 1.04, 1.11, 1.17, 1.23, 1.29, 1.36, 1.42, 1.48, 1.54, 1.87, 2.19,
441  2.51, 2.83, 3.15, 3.48, 3.80, 4.12, 3.88, 3.65, 3.41, 3.18, 2.94,
442  2.71, 2.47, 2.24, 2.16, 2.08, 2.00, 1.93, 1.85, 1.77, 1.69, 1.62,
443  1.54},
444 
445  { 0.70, 0.76, 0.81, 0.87, 0.92, 0.98, 1.03, 1.09, 1.14, 1.20, 1.56, 1.93,
446  2.29, 2.66, 3.02, 3.39, 3.75, 4.12, 3.93, 3.74, 3.54, 3.35, 3.16,
447  2.97, 2.78, 2.58, 2.50, 2.41, 2.33, 2.24, 2.16, 2.07, 1.99, 1.90,
448  1.82}
449  };
450 
451  int days[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
452 
453  // Set up for the time interpolation:
454 
455  if (year % 4 == 0)
456  days[1] = days[1] + 1;
457 
458  daymid = days[month] / 2.0;
459 
460  if (day >= daymid) {
461  m1 = month;
462  m2 = m1 + 1;
463  if (m2 > 11)
464  m2 = 0;
465  t = day;
466  } else {
467  m2 = month;
468  m1 = m2 - 1;
469  if (m1 < 0)
470  m1 = 11;
471  t = days[m1] + day;
472  }
473 
474  t1 = days[m1] / 2.0;
475  t2 = days[m1] + (days[m2] / 2.0);
476 
477  // Perform the lat. interpolation, and further set up for time interp.:
478 
479  if (lat >= 85.0) {
480  d1 = tabwv[m1][0];
481  d2 = tabwv[m2][0];
482  } else if (lat <= (-85.0)) {
483  d1 = tabwv[m1][34];
484  d2 = tabwv[m2][34];
485  }
486  else {
487  i1 = 16 - (int) (lat / 5.0);
488  i2 = i1 + 1;
489  fac = (tabwv[month][i2] - tabwv[month][i1]) / (-5.0);
490  difflat = lat - (90.0 - (i1 * 5.0));
491  d1 = tabwv[m1][i1] + (fac * difflat);
492  d2 = tabwv[m2][i1] + (fac * difflat);
493  }
494 
495  // Complete the time interpolation to get final estimate:
496 
497  fac2 = (d2 - d1) / (t2 - t1);
498 
499  watvap = d1 + (t - t1) * fac2;
500 
501  return watvap;
502 }
503 
504 /********************************************************************************
505  * Get Earth-Sun distance correction (VarSol)
506  *
507  * Calculation of the variability of the solar constant during the year.
508  * jday is the number of the day in the month
509  * dsol is a multiplicative factor to apply to the mean value of
510  * solar constant.
511  *
512  ********************************************************************************/
513 
514 float varsol(int32_t jday) {
515 
516  float om;
517  float dsol;
518 
519  om = (0.9856f * (float) (jday - 4)) * M_PI / 180.f;
520  dsol = 1.f / powf((1.f - .01673f * cosf(om)), 2.f);
521 
522  return dsol;
523 }
524 
525 /*************************************************************************************
526  * Compute the rise and set time for the given Julian day of year *
527  *************************************************************************************/
528 
529 void triseset(int32_t jday, float xlon, float xlat, float *trise, float *tset) {
530 
531  float fac, xlo, xla, xj, a1, a2, et, a3, delta, cosah, ah1, ah2, tsv1, tsv2,
532  tsm1, tsm2;
533 
534  fac = (float) M_PI / 180.f;
535 
536  xlo = xlon * fac;
537  xla = xlat * fac;
538  xj = (float) (jday - 1);
539 
540  a1 = (1.00554f * xj - 6.28306f) * fac;
541  a2 = (1.93946f * xj + 23.35089f) * fac;
542  et = -7.67825f * sinf(a1) - 10.09176f * sinf(a2);
543  a3 = (0.9683f * xj - 78.00878f) * fac;
544 
545  delta = 23.4856f * sinf(a3) * fac;
546  cosah = (-(sinf(xla) * sinf(delta))) / (cosf(xla) * cosf(delta));
547 
548  if ((cosah < -1.f) || (cosah > 1.f)) {
549  *trise = 0.0f;
550  *tset = 24.0f;
551  return;
552  }
553 
554  ah1 = acosf(cosah);
555  ah2 = -ah1;
556 
557  tsv1 = ah1 / (15.f * fac);
558  tsv2 = ah2 / (15.f * fac);
559  tsm1 = tsv1 + 12.f - et / 60.f;
560  tsm2 = tsv2 + 12.f - et / 60.f;
561 
562  *trise = tsm2 - (xlo / (15.f * fac));
563  *tset = tsm1 - (xlo / (15.f * fac));
564  return;
565 
566 }
567 
568 /********************************************************************************
569  * Get Sun position
570  ********************************************************************************/
571 
572 float get_solz(int32_t jday, float time, float lon, float lat) {
573  float asol;
574  double tsm, xla, xj, tet, a1, a2, a3, a4, a5, et, tsv, ah, b1, b2, b3, b4,
575  b5, b6, b7, delta, amuzero, elev, az, caz, azim, pi2;
576 
577  // solar position (zenithal angle asol,azimuthal angle phi0 in degrees)
578  // j is the day number in the year
579  //
580  // mean solar time (heure decimale) (aka decimal time)
581  double fac = M_PI / 180.;
582  tsm = time + lon / 15.;
583  xla = lat * fac;
584  xj = (float) (jday - 1);
585  tet = 2. * M_PI * xj / 365.;
586 
587  // time equation (in mn.dec)
588  a1 = 0.000075;
589  a2 = 0.001868;
590  a3 = 0.032077;
591  a4 = 0.014615;
592  a5 = 0.040849;
593  et = a1 + a2 * cos(tet) - a3 * sin(tet) - a4 * cos(2. * tet)
594  - a5 * sin(2. * tet);
595  et = et * 12. * 60. / M_PI;
596 
597  // true solar time
598 
599  tsv = tsm + et / 60.;
600  tsv -= 12.;
601 
602  // hour angle
603 
604  ah = tsv * 15. * fac;
605 
606  // solar declination (in radian)
607 
608  b1 = 0.006918;
609  b2 = 0.399912;
610  b3 = 0.070257;
611  b4 = 0.006758;
612  b5 = 0.000907;
613  b6 = 0.002697;
614  b7 = 0 - .001480;
615  delta = b1 - b2 * cos(tet) + b3 * sin(tet) - b4 * cos(2. * tet)
616  + b5 * sin(2. * tet) - b6 * cos(3. * tet) + b7 * sin(3.f * tet);
617 
618  // elevation,azimuth
619 
620  amuzero = sin(xla) * sin(delta) + cos(xla) * cos(delta) * cos(ah);
621  elev = asin(amuzero);
622  az = cos(delta) * sin(ah) / cos(elev);
623  if ((fabs(az) - 1.000) > 0.00000)
624  az = copysign(az, 1.);
625 
626  caz = (-cos(xla) * sin(delta) + sin(xla) * cos(delta) * cos(ah))
627  / cos(elev);
628  azim = asin(az);
629 
630  if (caz <= 0.)
631  azim = M_PI - azim;
632 
633  if (caz > 0. && az <= 0.)
634  azim = 2. * M_PI + azim;
635 
636  azim = azim + M_PI;
637  pi2 = 2.f * M_PI;
638 
639  if (azim > pi2)
640  azim -= pi2;
641 
642  elev /= fac;
643 
644  // conversion in degrees
645 
646  asol = (float) (90. - elev);
647  // phi0 = azim / fac;
648 
649  return asol;
650 
651 }
652 
653 /******************************************************************************
654  * Given inputs cos(SZ) and tauA500, interpolate the appropriate *
655  * surface albedo from the table. This is for the case where tau < 1.0 *
656  * Inputs are assumed to be in the valid ranges (the associated solar *
657  * zenith angle should be in range 0.0 - 87.134 degrees, and tau *
658  * should be in the range 0.0 - 0.99) *
659  ******************************************************************************/
660 
661 float interp_as_taulow(float csz, float tau) {
662 
663  int s1, s2, t1, t2, i;
664 
665  float stot, sdist, ttot, tdist, slope, alb1[2], as;
666 
667  float cszt[NS] = {0.05, 0.09, 0.15, 0.21, 0.27, 0.33, 0.39, 0.45, 0.52,
668  0.60, 0.68, 0.76, 0.84, 0.92, 1.00};
669 
670  float taut[NT] = {0.00, 0.05, 0.10, 0.16, 0.24, 0.35, 0.50, 0.70, 0.99};
671 
672  // Store the table here rather than read a new file. Fortran stores
673  // 2-d arrays internally by columns (col1, col2, col3 ...)
674 
675  float ast[NS][NT] = {
676  { 0.1592253, 0.1258933, 0.1088253, 0.0968594,
677  0.0880966, 0.0815492, 0.0766266, 0.0721790, 0.0679746},
678  {
679  0.1944972, 0.1637040, 0.1432702, 0.1240363, 0.1066057, 0.0922638,
680  0.0826043, 0.0757601, 0.0700167
681  },
682  { 0.1978650, 0.1785707,
683  0.1643212, 0.1482773, 0.1304572, 0.1119496, 0.0963537, 0.0839010,
684  0.0740948},
685  { 0.1785221, 0.1673987, 0.1589092, 0.1486409,
686  0.1359193, 0.1209508, 0.1061981, 0.0920763, 0.0793196},
687  {
688  0.1531512, 0.1473233, 0.1426633, 0.1366985, 0.1289257, 0.1188203,
689  0.1078114, 0.0957847, 0.0831087
690  },
691  { 0.1281988, 0.1255739,
692  0.1234628, 0.1204020, 0.1161201, 0.1101539, 0.1030820, 0.0943791,
693  0.0841410},
694  { 0.1061354, 0.1052812, 0.1048179, 0.1036617,
695  0.1016918, 0.0986726, 0.0948040, 0.0894514, 0.0822926},
696  {
697  0.0877530, 0.0881279, 0.0883850, 0.0884469, 0.0880587, 0.0870923,
698  0.0854952, 0.0826607, 0.0782380
699  },
700  { 0.0708031, 0.0720173,
701  0.0727886, 0.0736250, 0.0741367, 0.0746769, 0.0747173, 0.0741294,
702  0.0722523},
703  { 0.0567974, 0.0582123, 0.0593260, 0.0604172,
704  0.0615151, 0.0627740, 0.0639047, 0.0647193, 0.0650727},
705  {
706  0.0472026, 0.0485713, 0.0495830, 0.0507434, 0.0519943, 0.0536504,
707  0.0551176, 0.0566950, 0.0581553
708  },
709  { 0.0406631, 0.0419177,
710  0.0429259, 0.0440614, 0.0451823, 0.0467889, 0.0483670, 0.0501810,
711  0.0522433},
712  { 0.0366926, 0.0377500, 0.0384530, 0.0393431,
713  0.0404503, 0.0419569, 0.0434430, 0.0451987, 0.0473454},
714  {
715  0.0343793, 0.0352937, 0.0358116, 0.0364891, 0.0374246, 0.0385732,
716  0.0398924, 0.0414707, 0.0435983
717  },
718  { 0.0331561, 0.0337733,
719  0.0341567, 0.0346916, 0.0354239, 0.0364011, 0.0374280, 0.0387133,
720  0.0405543}
721  };
722 
723  // Set up interpolation:
724 
725  // Locate subcells in array to use:
726 
727  s1 = NS - 2;
728  s2 = NS - 1;
729  for (i = 0; i < NS - 1; i++) {
730  if (csz < cszt[i + 1]) {
731  s1 = i;
732  s2 = i + 1;
733  break;
734  }
735  }
736  t1 = NT - 2;
737  t2 = NT - 1;
738  for (i = 0; i < NT - 1; i++) {
739  if (tau < taut[i + 1]) {
740  t1 = i;
741  t2 = i + 1;
742  break;
743  }
744  }
745 
746  stot = cszt[s2] - cszt[s1];
747  sdist = csz - cszt[s1];
748  ttot = taut[t2] - taut[t1];
749  tdist = tau - taut[t1];
750 
751  slope = (ast[s2][t1] - ast[s1][t1]) / stot;
752  alb1[0] = ast[s1][t1] + (slope * sdist);
753  slope = (ast[s2][t2] - ast[s1][t2]) / stot;
754  alb1[1] = ast[s1][t2] + (slope * sdist);
755 
756  slope = (alb1[1] - alb1[0]) / ttot;
757  as = alb1[0] + (slope * tdist);
758 
759  return as;
760 }
761 
762 /********************************************************************************
763  * Given input cos(SZ), interpolate the appropriate surface albedo *
764  * from the table. This is for the case where tau > 1.0 *
765  * Input is assumed to be in the valid range (the associated solar *
766  * zenith angle should be in range 0.0 - 87.134 degrees) *
767  ********************************************************************************/
768 
769 float interp_as_tauhigh(float csz) {
770 
771  int s1, s2, i;
772 
773  float stot, sdist, slope, as;
774 
775  float cszt[NS] = {0.05, 0.09, 0.15, 0.21, 0.27, 0.33, 0.39, 0.45, 0.52,
776  0.60, 0.68, 0.76, 0.84, 0.92, 1.00};
777 
778  float ast[NS] = {0.0623166, 0.0630070, 0.0640739, 0.0650397, 0.0657760,
779  0.0661690, 0.0659415, 0.0651271, 0.0635101, 0.0610652, 0.0583470,
780  0.0556146, 0.0530646, 0.0509498, 0.0490149};
781 
782  // Set up interpolation:
783 
784  // Locate subcells in array to use:
785 
786  s1 = NS - 2;
787  s2 = NS - 1;
788  for (i = 0; i < NS; i++) {
789  if (csz < cszt[i + 1]) {
790  s1 = i;
791  s2 = i + 1;
792  break;
793  }
794  }
795 
796  stot = cszt[s2] - cszt[s1];
797  sdist = csz - cszt[s1];
798 
799  slope = (ast[s2] - ast[s1]) / stot;
800  as = ast[s1] + (slope * sdist);
801 
802  return as;
803 }
804 
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:240
data_t t[NROOTS+1]
Definition: decode_rs.h:77
int j
Definition: decode_rs.h:73
int32_t day
#define NT
Definition: par_utils.c:43
#define NMODEL
Definition: par_utils.c:41
#define NPHASE
Definition: par_utils.h:14
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
#define FALSE
Definition: rice.h:164
#define NULL
Definition: decode_rs.h:63
void read_aerosol_par(l2str *l2rec, int32_t nbands, float *tablewavelengths, float *tablephaseangles, float *tablealphas, float *tableomegas, float *tableaerphasefunc)
Definition: par_utils.c:174
#define TRUE
Definition: rice.h:165
float * lat
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start thereby resolving MODur00108 Changed header guard macro names to avoid reserved name resolving MODur00104 Maneuver list file for Terra satellite was updated to include data from Jecue DuChateu Maneuver list files for both Terra and Aqua were updated to include two maneuvers from recent IOT weekly reports The limits for Z component of angular momentum was and to set the fourth GEO scan quality flag for a scan depending upon whether or not it occurred during one of them Added _FillValue attributes to many and changed the fill value for the sector start times from resolving MODur00072 Writes boundingcoordinate metadata to L1A archived metadata For PERCENT *ECS change to treat as
Definition: HISTORY.txt:297
int32_t jday(int16_t i, int16_t j, int16_t k)
Definition: jday.c:4
void GetAerPhase(l2str *l2rec, int ip, int32_t nbands, float angstrom, float *aerphasefunc, float *omega, float *modelAngstrom)
Definition: par_utils.c:51
#define fac
double precision function f(R1)
Definition: tmd.lp.f:1454
float interp_as_tauhigh(float csz)
Definition: par_utils.c:769
data_t omega[NROOTS+1]
Definition: decode_rs.h:77
float32 slope[]
Definition: l2lists.h:30
const double delta
#define M_PI
Definition: pml_iop.h:15
instead the metadata field ProcessingEnvinronment is filled in from the output of a call to the POSIX compliant function uname from within the L1B code A small bug in L1B_Tables an incorrect comparison of RVS coefficients for TEBs to RVS coefficients for RSBs was being made This was replaced with a comparison between TEB coefficients This error never resulted in an incorrect RVS correction but did lead to recalculating the coefficients for each detector in a thermal band even if the coefficients were the same for all detectors To reduce to overall size of the reflective LUT HDF fill values were eliminated from all LUTs previously dimensioned where and where NUM_TIMES is the number of time dependent table pieces In Preprocess a small error where the trailing dropped scan counter was incremented when the leading dropped scan counter should have been was fixed This counter is internal only and is not yet used for any chiefly to casting of were added to make it LINUX compatible Output of code run on LINUX machines displays differences of at most scaled sector incalculable values of the Emissive calibration factor b1
Definition: HISTORY.txt:576
#define NS
Definition: par_utils.c:42
float xlon[LAC_PIXEL_NUM]
Definition: l1a_seawifs.c:93
void triseset(int32_t jday, float xlon, float xlat, float *trise, float *tset)
Definition: par_utils.c:529
float EstimateWatVap(int32_t year, int32_t month, int32_t day, float lat)
Definition: par_utils.c:381
float get_solz(int32_t jday, float time, float lon, float lat)
Definition: par_utils.c:572
int32_t nbands
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start time
Definition: HISTORY.txt:248
#define fabs(a)
Definition: misc.h:93
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
float EstimateDobson(int32_t year, int32_t month, int32_t day, float lat)
Definition: par_utils.c:250
float * lon
int i
Definition: decode_rs.h:71
float varsol(int32_t jday)
Definition: par_utils.c:514
float interp_as_taulow(float csz, float tau)
Definition: par_utils.c:661