A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
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 "par_utils.h"
38 
39 #include <math.h>
40 #include <stdio.h>
41 #include <stdlib.h>
42 
43 #define NMODEL 12
44 #define NS 15
45 #define NT 9
46 
47 static float *tablewavelengths;
48 static float *tablephaseangles;
49 static float *tablealphas; // angstrom exponent
50 static float *tableomegas; // single scattering albedo
51 static float *tableaerphasefunc;
52 
53 #define FILE_ERROR(str, name) { if (str.fid == -1) { \
54 printf("-Error--: could not open \"%s\" for reading.\n, See line %d in %s. Exiting...", name, __LINE__, __FILE__); \
55 exit(EXIT_FAILURE);\
56 }}
57 
58 inline float kasten_equation(float solz) {
59  float CosSZ = cos(solz * M_PI / 180.0);
60  return (1.003198 * CosSZ + 0.101632) /
61  (CosSZ * CosSZ + 0.090560 * CosSZ + 0.003198);
62 }
63 
64 size_t search(const float *arr, size_t s, size_t e, float val, size_t *i_s,
65  size_t *i_e) {
66  const bool acsending = arr[s] < arr[e];
67  if (acsending) {
68  if (val >= arr[e]) {
69  *i_s = e;
70  *i_e = e;
71  return e;
72  }
73  if (val <= arr[s]) {
74  *i_s = s;
75  *i_e = s;
76  return s;
77  }
78  } else {
79  if (val <= arr[e]) {
80  *i_s = e;
81  *i_e = e;
82  return e;
83  }
84  if (val >= arr[s]) {
85  *i_s = s;
86  *i_e = s;
87  return s;
88  }
89  }
90  while (e - s > 1) {
91  size_t m = (s + e) / 2; // compute median
92  if (acsending) {
93  if (arr[m] <= val) {
94  s = m;
95  } else {
96  e = m;
97  }
98  } else {
99  if (arr[m] <= val) {
100  e = m;
101  } else {
102  s = m;
103  }
104  }
105  }
106  {
107  *i_s = s;
108  *i_e = s + 1;
109  }
110  return s;
111 }
112 
113 size_t get_reduced_product(size_t n_dims, const size_t *dims) {
114  size_t result = 1;
115  for (size_t i = 0; i < n_dims; i++) {
116  result *= dims[i];
117  }
118  return result;
119 }
120 
121 size_t get_hyper_cube_index(size_t n_dims, const size_t *dims,
122  const size_t *indexes) {
123  if (n_dims > 1)
124  return indexes[0] * get_reduced_product(n_dims - 1, dims + 1) +
125  get_hyper_cube_index(n_dims - 1, dims + 1, indexes + 1);
126  else
127  return indexes[0];
128 }
129 
130 // 2^n_dims, float width, float diff
131 float hyper_cube_linear_interp(size_t n_dims, float *hyper_cube, float *width) {
132  if (n_dims == 1) // hyper_cube[2]
133  {
134  float ans = (hyper_cube[1] - hyper_cube[0]) * width[0] + hyper_cube[0];
135  return ans;
136  } else {
137  size_t total_index = pow(2, n_dims - 1);
138  float *small_hypercube = (float *) calloc(total_index, sizeof(float));
139  for (size_t ind = 0; ind < total_index; ind++) {
140  small_hypercube[ind] =
141  (hyper_cube[ind + total_index] - hyper_cube[ind]) * width[0] +
142  hyper_cube[ind];
143  // {
144  // std::cout << "small hypercube " << n_dims << " " << ind << "
145  // "
146  // << small_hypercube[ind] << " "
147  // << hyper_cube[ind + total_index] << " "
148  // << hyper_cube[ind] << " " << width[0] << std::endl;
149  // }
150  }
151  float ans =
152  hyper_cube_linear_interp(n_dims - 1, small_hypercube, width + 1);
153  free(small_hypercube);
154  return ans;
155  }
156 }
157 
158 float interpnd(size_t n_dims, const size_t *dims, const float *point,
159  float **grid, const float *lut) {
160  size_t *st_pts = (size_t *) calloc(n_dims, sizeof(size_t));
161  size_t *end_pts = (size_t *) calloc(n_dims, sizeof(size_t));;
162  for (size_t i_dim = 0; i_dim < n_dims; i_dim++) {
163  size_t dim_size = dims[i_dim];
164  size_t st = 0;
165  size_t end = dim_size - 1;
166  search(grid[i_dim], st, end, point[i_dim], &st, &end);
167  st_pts[i_dim] = st;
168  end_pts[i_dim] = end;
169  }
170  size_t total_index = pow(2, n_dims);
171  float *hypercube = (float *) calloc(total_index, sizeof(float));
172  float *width = (float *) calloc(n_dims, sizeof(float));
173  size_t *indexes = (size_t *) calloc(n_dims, sizeof(size_t));
174  for (size_t i_dim = 0; i_dim < n_dims; i_dim++) {
175  size_t st = st_pts[i_dim];
176  size_t end = end_pts[i_dim];
177  if (st != end)
178  width[i_dim] = (point[i_dim] - grid[i_dim][st]) /
179  (grid[i_dim][end] - grid[i_dim][st]);
180  else
181  width[i_dim] = 0.0;
182  // { std::cout << "width[i_dim]" << width[i_dim] << std::endl; }
183  }
184 
185  for (size_t i = 0; i < total_index; i++) {
186  size_t hyper_cube_index = 0;
187  for (size_t i_dim = 0; i_dim < n_dims; i_dim++) {
188  size_t bit = (i >> i_dim) & 1;
189  if (bit) {
190  indexes[i_dim] = end_pts[i_dim];
191  hyper_cube_index += pow(2, n_dims - i_dim - 1);
192  } else {
193  indexes[i_dim] = st_pts[i_dim];
194  }
195  // {
196  // std::cout << "index " << i_dim << " i_dim " << indexes[i_dim]
197  // << std::endl;
198  // }
199  }
200  size_t index_lut = get_hyper_cube_index(n_dims, dims, indexes);
201  hypercube[hyper_cube_index] = lut[index_lut];
202  // {
203  // std::cout << i << " " << hyper_cube_index << " " << index_lut <<
204  // " "
205  // << hypercube[hyper_cube_index] << std::endl;
206  // }
207  }
208  float ans = hyper_cube_linear_interp(n_dims, hypercube, width);
209  free(st_pts);
210  free(end_pts);
211  free(hypercube);
212  free(width);
213  free(indexes);
214  return ans;
215 }
216 
217 void get_nc_dim(const idDS *struct_nc, const char *name, size_t *len) {
218  int dim_id;
219  int status = nc_inq_dimid((*struct_nc).fid, name, &dim_id);
220  if (status == 0) status = nc_inq_dimlen((*struct_nc).fid, dim_id, len);
221  if (status != 0) {
222  printf("--Error--: Dimension %s could not be read. See line %d in %s. Exiting...", name, __LINE__, __FILE__);
223  exit(EXIT_FAILURE);
224  }
225 }
226 
227 void get_nc_var(const idDS *struct_nc, const char *name, float *var,
228  const size_t *start, const size_t *count) {
229  int var_id;
230  int status = nc_inq_varid((*struct_nc).fid, name, &var_id);
231  if (status == 0)
232  status = nc_get_vara_float((*struct_nc).fid, var_id, start, count, var);
233  if (status != 0) {
234  printf("--Error--: Variable %s could not be read. See line %d in %s. Exiting...", name, __LINE__, __FILE__);
235  exit(EXIT_FAILURE);
236  }
237 }
238 
244 void get_luts_data(l2str *l2rec, luts_par *luts_data) {
245  const char *dataroot = getenv("OCDATAROOT");
246  char lut_rho_path[FILENAME_MAX] = "";
247  char lut_tg_path[FILENAME_MAX] = "";
248  char lut_td_path[FILENAME_MAX] = "";
249  char lut_dobson_path[FILENAME_MAX] = "";
250  char lut_watvap_path[FILENAME_MAX] = "";
251  char lut_surfoceanalbedo[FILENAME_MAX] = "";
252  char luts_scalar_par_path[FILENAME_MAX] = "";
253  char luts_scalar_inst_par_path[FILENAME_MAX] = "";
254  sprintf(lut_rho_path, "%s/common/LUT_rho.nc", dataroot);
255  sprintf(lut_td_path, "%s/common/LUT_td.nc", dataroot);
256  sprintf(lut_tg_path, "%s/common/LUT_tg.nc", dataroot);
257  sprintf(lut_dobson_path, "%s/common/LUT_dobson.nc", dataroot);
258  sprintf(lut_watvap_path, "%s/common/LUT_watvap.nc", dataroot);
259  sprintf(lut_surfoceanalbedo, "%s/common/LUT_ocean_surface_albedo.nc",
260  dataroot);
261  sprintf(luts_scalar_par_path, "%s/common/LUT_scalar_PAR.nc", dataroot);
262  sprintf(luts_scalar_inst_par_path, "%s/common/LUT_scalar_inst_PAR.nc", dataroot);
263  int16_t year, month, mday, doy;
264  double sec;
265  unix2yds(l2rec->l1rec->scantime, &year, &doy, &sec);
266  unix2ymds(l2rec->l1rec->scantime, &year, &month, &mday, &sec);
267  // Reading rho lut
268  {
269  idDS luts_rho_struct = openDS(lut_rho_path);
270  FILE_ERROR(luts_rho_struct, lut_rho_path)
271  size_t solar_zenith_angle, view_zenith_angle, relative_azimuth_angle,
272  angstrom_coefficients, optical_depth_at_550_nm, wavelength;
273  get_nc_dim(&luts_rho_struct, "solar_zenith_angle", &solar_zenith_angle);
274  get_nc_dim(&luts_rho_struct, "view_zenith_angle", &view_zenith_angle);
275  get_nc_dim(&luts_rho_struct, "relative_azimuth_angle",
276  &relative_azimuth_angle);
277  get_nc_dim(&luts_rho_struct, "angstrom_coefficients",
278  &angstrom_coefficients);
279  get_nc_dim(&luts_rho_struct, "optical_depth_at_550_nm",
280  &optical_depth_at_550_nm);
281  get_nc_dim(&luts_rho_struct, "wavelength", &wavelength);
282  size_t start[] = {0, 0, 0, 0, 0, 0};
283  size_t count[] = {solar_zenith_angle, view_zenith_angle,
284  relative_azimuth_angle, angstrom_coefficients,
285  optical_depth_at_550_nm, wavelength};
286  const size_t size_all = solar_zenith_angle * view_zenith_angle *
287  relative_azimuth_angle * angstrom_coefficients *
288  optical_depth_at_550_nm * wavelength;
289  luts_data->lut_rho = calloc(size_all, sizeof(float));
290 
291  luts_data->rhodims.dimsolar_zenith_angle = solar_zenith_angle;
292  luts_data->rhodims.dimview_zenith_angle = view_zenith_angle;
293  luts_data->rhodims.dimrelative_azimuth_angle = relative_azimuth_angle;
294  luts_data->rhodims.dimangstrom_coefficients = angstrom_coefficients;
295  luts_data->rhodims.dimoptical_depth_at_550_nm = optical_depth_at_550_nm;
296  luts_data->rhodims.dimwavelength = wavelength;
297 
298  luts_data->rhodims.solar_zenith_angle =
299  calloc(solar_zenith_angle, sizeof(float));
300  luts_data->rhodims.view_zenith_angle =
301  calloc(view_zenith_angle, sizeof(float));
302  luts_data->rhodims.relative_azimuth_angle =
303  calloc(relative_azimuth_angle, sizeof(float));
304  luts_data->rhodims.angstrom_coefficients =
305  calloc(angstrom_coefficients, sizeof(float));
306  luts_data->rhodims.optical_depth_at_550_nm =
307  calloc(optical_depth_at_550_nm, sizeof(float));
308  luts_data->rhodims.wavelength = calloc(wavelength, sizeof(float));
309  get_nc_var(&luts_rho_struct, "lut_rho", luts_data->lut_rho, start,
310  count);
311  size_t start_dim[] = {0};
312  size_t count_dim[] = {0};
313  count_dim[0] = solar_zenith_angle;
314  get_nc_var(&luts_rho_struct, "solar_zenith_angle",
315  luts_data->rhodims.solar_zenith_angle, start_dim, count_dim);
316 
317  count_dim[0] = view_zenith_angle;
318  get_nc_var(&luts_rho_struct, "view_zenith_angle",
319  luts_data->rhodims.view_zenith_angle, start_dim, count_dim);
320 
321  count_dim[0] = relative_azimuth_angle;
322  get_nc_var(&luts_rho_struct, "relative_azimuth_angle",
323  luts_data->rhodims.relative_azimuth_angle, start_dim,
324  count_dim);
325 
326  count_dim[0] = angstrom_coefficients;
327  get_nc_var(&luts_rho_struct, "angstrom_coefficients",
328  luts_data->rhodims.angstrom_coefficients, start_dim,
329  count_dim);
330 
331  count_dim[0] = optical_depth_at_550_nm;
332  get_nc_var(&luts_rho_struct, "optical_depth_at_550_nm",
333  luts_data->rhodims.optical_depth_at_550_nm, start_dim,
334  count_dim);
335 
336  count_dim[0] = wavelength;
337  get_nc_var(&luts_rho_struct, "wavelength",
338  luts_data->rhodims.wavelength, start_dim, count_dim);
339 
340  endDS(luts_rho_struct);
341  }
342  // Reading tg lut
343  {
344  idDS luts_tg_struct = openDS(lut_tg_path);
345  FILE_ERROR(luts_tg_struct, lut_tg_path)
346  size_t wavelength, air_mass, ozone_concentration, water_vapor_pressure;
347  get_nc_dim(&luts_tg_struct, "wavelength", &wavelength);
348  get_nc_dim(&luts_tg_struct, "air_mass", &air_mass);
349  get_nc_dim(&luts_tg_struct, "ozone_concentration",
350  &ozone_concentration);
351  get_nc_dim(&luts_tg_struct, "water_vapor_pressure",
352  &water_vapor_pressure);
353  size_t start[] = {0, 0, 0, 0};
354  size_t count[] = {wavelength, air_mass, water_vapor_pressure,
355  ozone_concentration};
356  const size_t size_all =
357  wavelength * air_mass * ozone_concentration * water_vapor_pressure;
358  luts_data->tgdims.dimair_mass = air_mass;
359  luts_data->tgdims.dimozone_concentration = ozone_concentration;
360  luts_data->tgdims.dimwavelength = wavelength;
361  luts_data->tgdims.dimwater_vapor_pressure = water_vapor_pressure;
362  luts_data->lut_tg = calloc(size_all, sizeof(float));
363  luts_data->tgdims.air_mass = calloc(air_mass, sizeof(float));
364  luts_data->tgdims.wavelength = calloc(wavelength, sizeof(float));
365  luts_data->tgdims.ozone_concentration =
366  calloc(ozone_concentration, sizeof(float));
367  luts_data->tgdims.water_vapor_pressure =
368  calloc(water_vapor_pressure, sizeof(float));
369 
370  get_nc_var(&luts_tg_struct, "lut_tg", luts_data->lut_tg, start, count);
371 
372  size_t start_dim[] = {0};
373  size_t count_dim[] = {0};
374 
375  count_dim[0] = wavelength;
376  get_nc_var(&luts_tg_struct, "wavelength", luts_data->tgdims.wavelength,
377  start_dim, count_dim);
378 
379  count_dim[0] = water_vapor_pressure;
380  get_nc_var(&luts_tg_struct, "water_vapor_pressure",
381  luts_data->tgdims.water_vapor_pressure, start_dim,
382  count_dim);
383 
384  count_dim[0] = ozone_concentration;
385  get_nc_var(&luts_tg_struct, "ozone_concentration",
386  luts_data->tgdims.ozone_concentration, start_dim, count_dim);
387 
388  count_dim[0] = air_mass;
389  get_nc_var(&luts_tg_struct, "air_mass", luts_data->tgdims.air_mass,
390  start_dim, count_dim);
391  endDS(luts_tg_struct);
392  }
393  // Reading td lut
394  {
395  idDS luts_td_struct = openDS(lut_td_path);
396  FILE_ERROR(luts_td_struct, lut_td_path)
397  size_t wavelength, air_mass, angstrom_coefficients,
398  optical_depth_at_550_nm;
399  get_nc_dim(&luts_td_struct, "wavelength", &wavelength);
400  get_nc_dim(&luts_td_struct, "air_mass", &air_mass);
401  get_nc_dim(&luts_td_struct, "angstrom_coefficients",
402  &angstrom_coefficients);
403  get_nc_dim(&luts_td_struct, "optical_depth_at_550_nm",
404  &optical_depth_at_550_nm);
405 
406  size_t start[] = {0, 0, 0, 0};
407  size_t count[] = {wavelength, air_mass, angstrom_coefficients,
408  optical_depth_at_550_nm};
409  const size_t size_all = wavelength * air_mass * angstrom_coefficients *
410  optical_depth_at_550_nm;
411 
412  luts_data->tddims.dimair_mass = air_mass;
413  luts_data->tddims.dimoptical_depth_at_550_nm = optical_depth_at_550_nm;
414  luts_data->tddims.dimwavelength = wavelength;
415  luts_data->tddims.dimangstrom_coefficients = angstrom_coefficients;
416  luts_data->lut_td = calloc(size_all, sizeof(float));
417 
418  luts_data->tddims.air_mass = calloc(air_mass, sizeof(float));
419  luts_data->tddims.wavelength = calloc(wavelength, sizeof(float));
420  luts_data->tddims.angstrom_coefficients =
421  calloc(angstrom_coefficients, sizeof(float));
422  luts_data->tddims.optical_depth_at_550_nm =
423  calloc(optical_depth_at_550_nm, sizeof(float));
424  get_nc_var(&luts_td_struct, "lut_td", luts_data->lut_td, start, count);
425 
426  size_t start_dim[] = {0};
427  size_t count_dim[] = {0};
428 
429  count_dim[0] = wavelength;
430  get_nc_var(&luts_td_struct, "wavelength", luts_data->tddims.wavelength,
431  start_dim, count_dim);
432 
433  count_dim[0] = angstrom_coefficients;
434  get_nc_var(&luts_td_struct, "angstrom_coefficients",
435  luts_data->tddims.angstrom_coefficients, start_dim,
436  count_dim);
437 
438  count_dim[0] = optical_depth_at_550_nm;
439  get_nc_var(&luts_td_struct, "optical_depth_at_550_nm",
440  luts_data->tddims.optical_depth_at_550_nm, start_dim,
441  count_dim);
442 
443  count_dim[0] = air_mass;
444  get_nc_var(&luts_td_struct, "air_mass", luts_data->tddims.air_mass,
445  start_dim, count_dim);
446  endDS(luts_td_struct);
447  }
448 
449  // Reading Ozone Conc. Lut
450  {
451  idDS luts_dobson_struct = openDS(lut_dobson_path);
452  FILE_ERROR(luts_dobson_struct, lut_dobson_path)
453  size_t days, latitude;
454  get_nc_dim(&luts_dobson_struct, "days", &days);
455  get_nc_dim(&luts_dobson_struct, "latitude", &latitude);
456  luts_data->ozonedims.dimlatitude = latitude;
457  luts_data->ozonedims.dimdays = days;
458  size_t start[] = {
459  0,
460  0,
461  };
462  size_t count[] = {days, latitude};
463  const size_t size_all = days * latitude;
464  luts_data->lut_dobson = calloc(size_all, sizeof(float));
465  luts_data->ozonedims.latitude = calloc(latitude, sizeof(float));
466  luts_data->ozonedims.days = calloc(days, sizeof(float));
467  get_nc_var(&luts_dobson_struct, "Dobson", luts_data->lut_dobson, start,
468  count);
469  size_t start_dim[] = {0};
470  size_t count_dim[] = {0};
471  count_dim[0] = latitude;
472  get_nc_var(&luts_dobson_struct, "latitude",
473  luts_data->ozonedims.latitude, start_dim, count_dim);
474  count_dim[0] = days;
475  if (year % 4 == 0)
476  get_nc_var(&luts_dobson_struct, "days_leap",
477  luts_data->ozonedims.days, start_dim, count_dim);
478  else
479  get_nc_var(&luts_dobson_struct, "days_no_leap",
480  luts_data->ozonedims.days, start_dim, count_dim);
481  endDS(luts_dobson_struct);
482  }
483  // Reading Wat. Vapor. Lut.
484  {
485  idDS luts_watvap_struct = openDS(lut_watvap_path);
486  FILE_ERROR(luts_watvap_struct, lut_watvap_path)
487  size_t days, latitude;
488  get_nc_dim(&luts_watvap_struct, "days", &days);
489  get_nc_dim(&luts_watvap_struct, "latitude", &latitude);
490  luts_data->watvapdims.dimlatitude = latitude;
491  luts_data->watvapdims.dimdays = days;
492  size_t start[] = {
493  0,
494  0,
495  };
496  size_t count[] = {days, latitude};
497  const size_t size_all = days * latitude;
498  luts_data->lut_watvap = calloc(size_all, sizeof(float));
499  luts_data->watvapdims.latitude = calloc(latitude, sizeof(float));
500  luts_data->watvapdims.days = calloc(days, sizeof(float));
501  get_nc_var(&luts_watvap_struct, "WatVap", luts_data->lut_watvap, start,
502  count);
503  size_t start_dim[] = {0};
504  size_t count_dim[] = {0};
505  count_dim[0] = latitude;
506  get_nc_var(&luts_watvap_struct, "latitude",
507  luts_data->watvapdims.latitude, start_dim, count_dim);
508  count_dim[0] = days;
509  if (year % 4 == 0)
510  get_nc_var(&luts_watvap_struct, "days_leap",
511  luts_data->watvapdims.days, start_dim, count_dim);
512  else
513  get_nc_var(&luts_watvap_struct, "days_no_leap",
514  luts_data->watvapdims.days, start_dim, count_dim);
515  endDS(luts_watvap_struct);
516  }
517  // Reading Surface Ocean Albedo Lut
518  {
519  idDS luts_soa_struct = openDS(lut_surfoceanalbedo);
520  FILE_ERROR(luts_soa_struct, lut_surfoceanalbedo)
521  size_t dim_jwvl, dim_wv1, dim_wv3, dim_wlt_reft;
522  get_nc_dim(&luts_soa_struct, "dim_jwvl", &dim_jwvl);
523  get_nc_dim(&luts_soa_struct, "dim_wv1", &dim_wv1);
524  get_nc_dim(&luts_soa_struct, "dim_wv3", &dim_wv3);
525  get_nc_dim(&luts_soa_struct, "dim_wlt_reft", &dim_wlt_reft);
526 
527  luts_data->soa_lut.dim_jwvl = dim_jwvl;
528  luts_data->soa_lut.dim_wv1 = dim_wv1;
529  luts_data->soa_lut.dim_wv3 = dim_wv3;
530  luts_data->soa_lut.dim_wlt_reft = dim_wlt_reft;
531 
532  luts_data->soa_lut.jwl = (float *) calloc(dim_jwvl, sizeof(float));
533  luts_data->soa_lut.achl = (float *) calloc(dim_jwvl, sizeof(float));
534 
535  luts_data->soa_lut.wv1 = (float *) calloc(dim_wv1, sizeof(float));
536  luts_data->soa_lut.bw1 = (float *) calloc(dim_wv1, sizeof(float));
537 
538  luts_data->soa_lut.wv3 = (float *) calloc(dim_wv3, sizeof(float));
539  luts_data->soa_lut.aw3 = (float *) calloc(dim_wv3, sizeof(float));
540 
541  luts_data->soa_lut.wlt_reft =
542  (float *) calloc(dim_wlt_reft, sizeof(float));
543  luts_data->soa_lut.reft = (float *) calloc(dim_wlt_reft, sizeof(float));
544 
545  size_t start_dim[] = {0};
546  size_t count_dim[] = {0};
547 
548  count_dim[0] = dim_jwvl;
549  get_nc_var(&luts_soa_struct, "jwvl", luts_data->soa_lut.jwl, start_dim,
550  count_dim);
551  get_nc_var(&luts_soa_struct, "achl", luts_data->soa_lut.achl, start_dim,
552  count_dim);
553 
554  count_dim[0] = dim_wv1;
555  get_nc_var(&luts_soa_struct, "wv1", luts_data->soa_lut.wv1, start_dim,
556  count_dim);
557  get_nc_var(&luts_soa_struct, "bw1", luts_data->soa_lut.bw1, start_dim,
558  count_dim);
559 
560  count_dim[0] = dim_wv3;
561  get_nc_var(&luts_soa_struct, "wv3", luts_data->soa_lut.wv3, start_dim,
562  count_dim);
563  get_nc_var(&luts_soa_struct, "aw3", luts_data->soa_lut.aw3, start_dim,
564  count_dim);
565 
566  count_dim[0] = dim_wlt_reft;
567  get_nc_var(&luts_soa_struct, "wlt_reft", luts_data->soa_lut.wlt_reft,
568  start_dim, count_dim);
569  get_nc_var(&luts_soa_struct, "reft", luts_data->soa_lut.reft, start_dim,
570  count_dim);
571  endDS(luts_soa_struct);
572  }
573  // Reading Scalar Par and Mean Cosine Luts
574  {
575  idDS luts_scalar_par_struct = openDS(luts_scalar_par_path);
576  FILE_ERROR(luts_scalar_par_struct, luts_scalar_par_path)
577  size_t dim_wind_speed, dim_doy, dim_latitude;
578  get_nc_dim(&luts_scalar_par_struct, "wind_speed", &dim_wind_speed);
579  get_nc_dim(&luts_scalar_par_struct, "doy", &dim_doy);
580  get_nc_dim(&luts_scalar_par_struct, "latitude", &dim_latitude);
581  luts_data->scalar_luts.dim_latitude = dim_latitude;
582  luts_data->scalar_luts.dim_wind_speed = dim_wind_speed;
583  luts_data->scalar_luts.dim_doy = dim_doy;
584  luts_data->scalar_luts.latitude =
585  (float *) calloc(dim_latitude, sizeof(float));
586  luts_data->scalar_luts.doy = (float *) calloc(dim_doy, sizeof(float));
587  luts_data->scalar_luts.wind_speed =
588  (float *) calloc(dim_wind_speed, sizeof(float));
589  size_t start_dim[] = {0};
590  size_t count_dim[] = {0};
591  count_dim[0] = dim_latitude;
592  get_nc_var(&luts_scalar_par_struct, "latitude",
593  luts_data->scalar_luts.latitude, start_dim, count_dim);
594  count_dim[0] = dim_wind_speed;
595  get_nc_var(&luts_scalar_par_struct, "wind_speed",
596  luts_data->scalar_luts.wind_speed, start_dim, count_dim);
597  count_dim[0] = dim_doy;
598  get_nc_var(&luts_scalar_par_struct, "doy", luts_data->scalar_luts.doy,
599  start_dim, count_dim);
600 
601  size_t start[] = {0, 0, 0};
602  size_t count[] = {dim_wind_speed, dim_doy, dim_latitude};
603  const size_t size_all = dim_wind_speed * dim_doy * dim_latitude;
604  luts_data->scalar_luts.CF = (float *) calloc(size_all, sizeof(float));
605  luts_data->scalar_luts.mu_cosine =
606  (float *) calloc(size_all, sizeof(float));
607  luts_data->scalar_luts.mu_cosine_c =
608  (float *) calloc(size_all, sizeof(float));
609  luts_data->scalar_luts.PARo = (float *) calloc(size_all, sizeof(float));
610  luts_data->scalar_luts.PARo_c =
611  (float *) calloc(size_all, sizeof(float));
612  luts_data->scalar_luts.PARc_above =
613  (float *) calloc(size_all, sizeof(float));
614  get_nc_var(&luts_scalar_par_struct, "CF", luts_data->scalar_luts.CF,
615  start, count);
616  get_nc_var(&luts_scalar_par_struct, "mu_cosine",
617  luts_data->scalar_luts.mu_cosine, start, count);
618  get_nc_var(&luts_scalar_par_struct, "mu_cosine_c",
619  luts_data->scalar_luts.mu_cosine_c, start, count);
620  get_nc_var(&luts_scalar_par_struct, "PARo", luts_data->scalar_luts.PARo,
621  start, count);
622  get_nc_var(&luts_scalar_par_struct, "PARo_c",
623  luts_data->scalar_luts.PARo_c, start, count);
624  get_nc_var(&luts_scalar_par_struct, "PARc_above",
625  luts_data->scalar_luts.PARc_above, start, count);
626 
627  }
628  // Reading Scalar Inst Par LUT
629  {
630  idDS luts_scalar_par_inst_struct = openDS(luts_scalar_inst_par_path);
631  FILE_ERROR(luts_scalar_par_inst_struct, luts_scalar_inst_par_path)
632  size_t dim_wind_speed, dim_solz, dim_cot;
633  get_nc_dim(&luts_scalar_par_inst_struct, "wind_speed", &dim_wind_speed);
634  get_nc_dim(&luts_scalar_par_inst_struct, "solz", &dim_solz);
635  get_nc_dim(&luts_scalar_par_inst_struct, "cot", &dim_cot); // cot is for cf grid
636  luts_data->scalar_inst_luts.dim_wind_speed = dim_wind_speed;
637  luts_data->scalar_inst_luts.dim_solz = dim_solz;
638  luts_data->scalar_inst_luts.dim_cot = dim_cot;
639  luts_data->scalar_inst_luts.wind_speed = (float *) calloc(dim_wind_speed, sizeof(float));
640  luts_data->scalar_inst_luts.cos_solz = (float *) calloc(dim_solz, sizeof(float));
641  luts_data->scalar_inst_luts.cot = (float *) calloc(dim_cot, sizeof(float));
642  luts_data->scalar_inst_luts.pard_m_cs = (float *) calloc(dim_wind_speed * dim_solz, sizeof(float));
643  luts_data->scalar_inst_luts.pard_p_cs = (float *) calloc(dim_wind_speed * dim_solz, sizeof(float));
644  luts_data->scalar_inst_luts.pard_m_oc = (float *) calloc(dim_wind_speed * dim_solz, sizeof(float));
645  luts_data->scalar_inst_luts.cf_pard_m = (float *) calloc(dim_wind_speed * dim_solz * dim_cot, sizeof(float));
646  luts_data->scalar_inst_luts.cf_pard_p = (float *) calloc(dim_wind_speed * dim_solz * dim_cot, sizeof(float));
647 
648  // reading 1d vars
649  {
650  size_t start[] = {0};
651  size_t count[] = {dim_wind_speed};
652  get_nc_var(&luts_scalar_par_inst_struct, "wind_speed", luts_data->scalar_inst_luts.wind_speed,
653  start, count);
654  count[0] = dim_solz;
655  get_nc_var(&luts_scalar_par_inst_struct, "cos_solz", luts_data->scalar_inst_luts.cos_solz,
656  start, count);
657  count[0] = dim_cot;
658  get_nc_var(&luts_scalar_par_inst_struct, "cot", luts_data->scalar_inst_luts.cot,
659  start, count);
660  }
661  // reading 2d and 3d vars
662  {
663  size_t start[] = {0, 0, 0};
664  size_t count[] = {dim_wind_speed, dim_solz, dim_cot};
665  get_nc_var(&luts_scalar_par_inst_struct, "PARo_below_clear_sky", luts_data->scalar_inst_luts.pard_m_cs,
666  start, count);
667  get_nc_var(&luts_scalar_par_inst_struct, "PARd_above_clear_sky", luts_data->scalar_inst_luts.pard_p_cs,
668  start, count);
669  get_nc_var(&luts_scalar_par_inst_struct, "PARo_below_overcast", luts_data->scalar_inst_luts.pard_m_oc,
670  start, count);
671 
672  get_nc_var(&luts_scalar_par_inst_struct, "CF_PARo_below", luts_data->scalar_inst_luts.cf_pard_m,
673  start, count);
674  get_nc_var(&luts_scalar_par_inst_struct, "CF_PARd_above", luts_data->scalar_inst_luts.cf_pard_p,
675  start, count);
676  }
677  }
678 }
679 
680 void GetAerPhase(l2str *l2rec, int ip, int32_t nbands, float angstrom,
681  float *aerphasefunc, float *omega, float *modelAngstrom) {
682  float phaseangle = l2rec->l1rec->scattang[ip];
683  float distx, dista, p1, p2;
684  int i, j, ix1 = 0, ix2, ia1 = 0, ia2, widx;
685 
686  static int firstCall = TRUE;
687 
688  // The data table is not in monotonically-increasing order, so here are the
689  // indices in proper order:
690 
691  int ii[12] = {1, 0, 5, 4, 8, 7, 3, 2, 6, 11, 10, 9};
692 
693  // Read the aerosol data file:
694 
695  if (firstCall) {
696  firstCall = FALSE;
697  tablewavelengths =
698  (float *) allocateMemory(nbands * sizeof(float), "tablewavelengths");
699  tablephaseangles =
700  (float *) allocateMemory(NPHASE * sizeof(float), "tablephaseangles");
701  tablealphas = (float *) allocateMemory(nbands * NMODEL * sizeof(float),
702  "tablealphas");
703  tableomegas = (float *) allocateMemory(nbands * NMODEL * sizeof(float),
704  "tableomegas");
705  tableaerphasefunc = (float *) allocateMemory(
706  NPHASE * nbands * NMODEL * sizeof(float), "tableaerphasefunc");
707  read_aerosol_par(l2rec, nbands, tablewavelengths, tablephaseangles,
708  tablealphas, tableomegas, tableaerphasefunc);
709  }
710 
711  // Find out which tablephaseangles-values in the table will be used:
712 
713  if (phaseangle < tablephaseangles[0]) {
714  // printf("Input phase angle too small, setting to 0. Was %f",
715  // phaseangle);
716  phaseangle = tablephaseangles[0];
717  ix1 = 0;
718  } else if (phaseangle >= tablephaseangles[NPHASE - 1]) {
719  // printf("Input phase angle too large, setting to 180. Was %f",
720  // phaseangle);
721  phaseangle = tablephaseangles[NPHASE - 1];
722  ix1 = NPHASE - 2;
723  } else {
724  for (i = NPHASE - 2; i > 0; i--) {
725  if (phaseangle >= tablephaseangles[i]) {
726  ix1 = i;
727  break;
728  }
729  }
730  }
731  ix2 = ix1 + 1;
732 
733  distx = (tablephaseangles[ix2] - phaseangle) /
734  (tablephaseangles[ix2] - tablephaseangles[ix1]);
735 
736  // Find out which models in the table to use:
737 
738  widx = windex(670.0, tablewavelengths, nbands);
739 
740  if (angstrom < tablealphas[widx * NMODEL + ii[0]]) {
741  angstrom = tablealphas[widx * NMODEL + ii[0]];
742  ia1 = 0;
743  } else if (angstrom >= tablealphas[widx * NMODEL + ii[NMODEL - 1]]) {
744  angstrom = tablealphas[widx * NMODEL + ii[NMODEL - 1]];
745  ia1 = NMODEL - 2;
746  } else {
747  for (j = NMODEL - 2; j > 0; j--) {
748  if (angstrom >= tablealphas[widx * NMODEL + ii[j]]) {
749  ia1 = j;
750  break;
751  }
752  }
753  }
754  ia2 = ia1 + 1;
755  dista = (tablealphas[widx * NMODEL + ii[ia2]] - angstrom) /
756  (tablealphas[widx * NMODEL + ii[ia2]] -
757  tablealphas[widx * NMODEL + ii[ia1]]);
758 
759  // Loop through the wavelengths and perform interpolations:
760 
761  for (i = 0; i < nbands; i++) {
762  omega[i] = (tableomegas[i * NMODEL + ii[ia1]] * dista) +
763  (tableomegas[i * NMODEL + ii[ia2]] * (1.0 - dista));
764 
765  modelAngstrom[i] = (tablealphas[i * NMODEL + ii[ia1]] * dista) +
766  (tablealphas[i * NMODEL + ii[ia2]] * (1.0 - dista));
767 
768  p1 = (tableaerphasefunc[i * NMODEL * NPHASE + ii[ia1] * NPHASE + ix1] *
769  distx) +
770  (tableaerphasefunc[i * NMODEL * NPHASE + ii[ia1] * NPHASE + ix2] *
771  (1.0 - distx));
772 
773  p2 = (tableaerphasefunc[i * NMODEL * NPHASE + ii[ia2] * NPHASE + ix1] *
774  distx) +
775  (tableaerphasefunc[i * NMODEL * NPHASE + ii[ia2] * NPHASE + ix2] *
776  (1.0 - distx));
777 
778  aerphasefunc[i] = (p1 * dista) + (p2 * (1.0 - dista));
779  }
780 
781  return;
782 }
783 
784 /********************************************************************************
785  *
786  * Description:
787  * This is to read the SeaWiFS aerosol data file for the PAR.
788  * Menghua Wang 6/22/99
789  * Parameters:
790  * angl, R, ---- scattering angles.
791  * omega0, R, --- single scattering albedo.
792  * alpha, R, --- Angstrom coefficient alpha(lambda,865).
793  * s11, R, --- aerosol phase function.
794  *
795  * SeaWiFS aerosol models:
796  * 1-2: Oceanic RH=90 and 99%
797  * 3-6: Maritime RH=50%, 70, 90, and 99%
798  * 7-9: Coastal RH=50, 90, and 99%
799  * 10-12: Tropospheric RH=50, 90, and 99%
800  *
801  *******************************************************************************/
802 void read_aerosol_par(l2str *l2rec, int32_t nbands, float *tablewavelengths,
803  float *tablephaseangles, float *tablealphas,
804  float *tableomegas, float *tableaerphasefunc) {
805  int imodel[NMODEL];
806 
807  char aerosolfilename[FILENAME_MAX] = "";
808  char *dataroot;
809  FILE *aerosoldata;
810  // Get OCDATAROOT path from user:
811  if ((dataroot = getenv("OCDATAROOT")) == NULL) {
812  printf("OCDATAROOT environment variable is not defined.\n");
813  return;
814  }
815  sprintf(aerosolfilename, "%s/%s/%s_aerosol_par.dat", dataroot,
816  sensorId2SensorDir(l2rec->l1rec->l1file->sensorID),
817  sensorId2SensorDir(l2rec->l1rec->l1file->sensorID));
818 
819  // Opens the aerosol data file for reading.
820  if ((aerosoldata = fopen(aerosolfilename, "r")) == NULL) {
821  printf("Error reading aerosol table for PAR from %s.\n",
822  aerosolfilename);
823  exit(EXIT_FAILURE);
824  }
825  printf("Loading aerosol properties for PAR from %s.\n", aerosolfilename);
826  int j, jj, iph;
827  int num_model;
828  for (jj = 0; jj < NPHASE; jj++) {
829  if ((fscanf(aerosoldata, "%f", &tablephaseangles[jj])) == 0) {
830  printf("Problem reading phase angles from %s.\n", aerosolfilename);
831  exit(EXIT_FAILURE);
832  }
833  }
834  for (j = 0; j < nbands; j++) {
835  for (iph = 0; iph < NMODEL; iph++) {
836  if ((fscanf(aerosoldata, "%d %d %f", &imodel[iph], &num_model,
837  &tablewavelengths[j])) == 0) {
838  printf("Problem reading model number and wavelength from %s.\n",
839  aerosolfilename);
840  exit(EXIT_FAILURE);
841  }
842  if ((fscanf(aerosoldata, "%f %f", &tableomegas[j * NMODEL + iph],
843  &tablealphas[j * NMODEL + iph])) == 0) {
844  printf("Problem reading SSA and Angstrom from %s.\n",
845  aerosolfilename);
846  exit(EXIT_FAILURE);
847  }
848  for (jj = 0; jj < NPHASE; jj++) {
849  if ((fscanf(aerosoldata, "%f",
850  &tableaerphasefunc[j * NMODEL * NPHASE +
851  iph * NPHASE + jj])) == 0) {
852  printf("Problem reading phase function values from %s.\n",
853  aerosolfilename);
854  exit(EXIT_FAILURE);
855  }
856  }
857  }
858  }
859  fclose(aerosoldata);
860 }
861 
862 /********************************************************************************
863  *
864  * Estimate Dobson units from climatology, given the month and
865  * latitude. Table has 12 columns, one per month, and 35 rows,
866  * from 85 N to -85, in steps of 5 deg lat.
867  * Perform 2d interpolation (latitude,month_day) to get estimate.
868  *
869  * Author: John McPherson <JMcPherson@ucsd.edu>
870  *
871  * Dobson look-up table from Keating et al., Ozone reference models
872  * for the middle atmosphere, Handbook for MAP, Vol. 31, 1-36, 1989
873  *
874  */
875 
876 float EstimateDobson(int32_t year, int32_t month, int32_t day, float lat) {
877  int i1, i2, m1, m2;
878  float daymid, d1, d2;
879  float t, t1, t2, fac, fac2, difflat;
880  float dobson;
881 
882  int tabdobson[12][35] = {
883  {395, 395, 395, 395, 395, 392, 390, 387, 376, 354, 322, 292,
884  269, 254, 248, 246, 247, 251, 255, 260, 266, 271, 277, 286,
885  295, 306, 319, 334, 344, 344, 338, 331, 324, 320, 316},
886 
887  {433, 433, 433, 436, 432, 428, 426, 418, 402, 374, 338, 303,
888  278, 261, 251, 246, 248, 250, 254, 258, 262, 265, 270, 278,
889  286, 294, 303, 313, 322, 325, 324, 317, 306, 299, 294},
890 
891  {467, 470, 460, 459, 451, 441, 433, 420, 401, 377, 347, 316,
892  291, 271, 260, 254, 254, 255, 257, 259, 261, 264, 269, 277,
893  284, 289, 296, 305, 312, 315, 317, 312, 305, 299, 295},
894 
895  {467, 465, 462, 455, 444, 431, 421, 410, 395, 373, 348, 325,
896  304, 287, 275, 267, 261, 259, 258, 259, 260, 263, 271, 278,
897  284, 289, 297, 306, 314, 318, 319, 313, 302, 302, 302},
898 
899  {411, 414, 416, 415, 410, 406, 402, 394, 382, 363, 342, 324,
900  307, 291, 279, 271, 264, 260, 258, 257, 258, 264, 271, 281,
901  291, 303, 312, 318, 322, 323, 322, 322, 322, 322, 322},
902 
903  {371, 371, 370, 368, 367, 372, 375, 372, 360, 341, 323, 311,
904  301, 290, 282, 275, 268, 263, 259, 256, 258, 264, 273, 289,
905  306, 319, 327, 328, 328, 337, 337, 337, 337, 337, 337},
906 
907  {333, 332, 332, 334, 338, 346, 350, 346, 335, 321, 310, 302,
908  296, 289, 284, 280, 274, 268, 262, 259, 261, 268, 279, 295,
909  315, 331, 340, 342, 338, 344, 340, 340, 340, 340, 340},
910 
911  {311, 308, 308, 313, 320, 327, 330, 326, 319, 310, 303, 298,
912  291, 286, 283, 281, 277, 273, 268, 264, 266, 274, 288, 306,
913  327, 343, 353, 355, 351, 339, 325, 307, 294, 294, 294},
914 
915  {283, 291, 302, 308, 312, 317, 318, 313, 307, 300, 295, 290,
916  284, 279, 279, 279, 278, 276, 272, 270, 273, 282, 295, 313,
917  333, 348, 360, 367, 368, 353, 324, 291, 267, 253, 230},
918 
919  {299, 299, 299, 309, 315, 317, 317, 312, 302, 291, 283, 280,
920  275, 270, 268, 267, 263, 263, 265, 269, 277, 287, 301, 317,
921  336, 354, 371, 387, 402, 402, 374, 333, 294, 274, 259},
922 
923  {314, 314, 314, 314, 332, 332, 327, 322, 311, 297, 284, 276,
924  270, 263, 261, 260, 258, 259, 264, 270, 278, 286, 298, 311,
925  323, 335, 350, 366, 381, 390, 388, 376, 357, 346, 341},
926 
927  {358, 358, 358, 358, 358, 358, 353, 349, 338, 320, 299, 281,
928  267, 256, 252, 251, 251, 253, 257, 264, 272, 279, 287, 297,
929  307, 318, 332, 347, 358, 365, 366, 364, 358, 356, 353}};
930 
931  int days[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
932 
933  // Set up for the time interpolation:
934 
935  if (year % 4 == 0) days[1] = days[1] + 1;
936 
937  daymid = days[month - 1] / 2.0;
938 
939  if (day >= daymid) {
940  m1 = month - 1;
941  m2 = m1 + 1;
942  if (m2 > 12) m2 = 1;
943  t = day;
944  } else {
945  m2 = month - 1;
946  m1 = m2 - 1;
947  if (m1 < 1) m1 = 12;
948  t = days[m1] + day;
949  }
950 
951  t1 = days[m1] / 2.0;
952  t2 = days[m1] + (days[m2] / 2.0);
953 
954  // Perform the lat. interpolation, and further set up for time interp.:
955 
956  if (lat >= 85.0) {
957  d1 = tabdobson[m1][1];
958  d2 = tabdobson[m2][1];
959  } else if (lat <= -85.0) {
960  d1 = tabdobson[m1][34];
961  d2 = tabdobson[m2][34];
962  } else {
963  i1 = 17 - (int) (lat / 5.0);
964  i2 = i1 + 1;
965  fac = (tabdobson[month][i2] - tabdobson[month][i1]) / (-5.0);
966  difflat = lat - (90.0 - (i1 * 5.0));
967  d1 = tabdobson[m1][i1] + (fac * difflat);
968  d2 = tabdobson[m2][i1] + (fac * difflat);
969  }
970 
971  // Complete the time interpolation to get final estimate:
972 
973  fac2 = (d2 - d1) / (t2 - t1);
974 
975  dobson = d1 + (t - t1) * fac2;
976  return dobson;
977 }
978 
979 /* ****************************************************************************
980  *
981  * Interpolate water vapor table, as for dobson (ozone).
982  * Table has 12 columns, one per month, and 35 rows,
983  * from 85 N to 85 S, in steps of 5 deg lat.
984  * winter peak (Jan 15), summer peak (July 15), north hemis
985  * Perform 2d interpolation (latitude,month_day) to get estimate.
986  *
987  * Tropical: if ( ABS( Lat ) .LE. 20.0 ) WatVap = 4.12
988  *
989  * Mid-latitude case: if ( ABS( Lat ) .LE. 60.0 ) Then
990  * if ( (North.And.Summer) .Or. (South.And.Winter) ) WatVap = 2.93
991  * else WatVap = 0.85
992  *
993  * Sub-arctic case:
994  * if ( (North.And.Summer) .Or. (South.And.Winter) ) WatVap = 2.10
995  * else WatVap = 0.42
996  *
997  * after McClatchey, R.A., et al, _Optical Properties of the Atmosphere_,
998  * AFCRL 71-0279 354, 91 pp., 1971
999  *
1000  */
1001 
1002 float EstimateWatVap(int32_t year, int32_t month, int32_t day, float lat) {
1003  int i1, i2, m1, m2;
1004  month -= 1;
1005  float daymid, d1, d2;
1006  float t, t1, t2, fac, fac2, difflat;
1007  float watvap;
1008 
1009  float tabwv[12][35] = {
1010  {0.42, 0.47, 0.52, 0.56, 0.61, 0.66, 0.71, 0.75, 0.80, 0.85, 1.26, 1.67,
1011  2.08, 2.48, 2.89, 3.30, 3.71, 4.12, 3.97, 3.82, 3.67, 3.53, 3.38, 3.23,
1012  3.08, 2.93, 2.84, 2.75, 2.65, 2.56, 2.47, 2.38, 2.28, 2.19, 2.10},
1013 
1014  {0.70, 0.76, 0.81, 0.87, 0.92, 0.98, 1.03, 1.09, 1.14, 1.20, 1.56, 1.93,
1015  2.29, 2.66, 3.02, 3.39, 3.75, 4.12, 3.93, 3.74, 3.54, 3.35, 3.16, 2.97,
1016  2.78, 2.58, 2.50, 2.41, 2.33, 2.24, 2.16, 2.07, 1.99, 1.90, 1.82},
1017 
1018  {0.98, 1.04, 1.11, 1.17, 1.23, 1.29, 1.36, 1.42, 1.48, 1.54, 1.87, 2.19,
1019  2.51, 2.83, 3.15, 3.48, 3.80, 4.12, 3.88, 3.65, 3.41, 3.18, 2.94, 2.71,
1020  2.47, 2.24, 2.16, 2.08, 2.00, 1.93, 1.85, 1.77, 1.69, 1.62, 1.54},
1021 
1022  {1.26, 1.33, 1.40, 1.47, 1.54, 1.61, 1.68, 1.75, 1.82, 1.89, 2.17, 2.45,
1023  2.73, 3.01, 3.28, 3.56, 3.84, 4.12, 3.84, 3.56, 3.28, 3.01, 2.73, 2.45,
1024  2.17, 1.89, 1.82, 1.75, 1.68, 1.61, 1.54, 1.47, 1.40, 1.33, 1.26},
1025 
1026  {1.54, 1.62, 1.69, 1.77, 1.85, 1.93, 2.00, 2.08, 2.16, 2.24, 2.47, 2.71,
1027  2.94, 3.18, 3.41, 3.65, 3.88, 4.12, 3.80, 3.48, 3.15, 2.83, 2.51, 2.19,
1028  1.87, 1.54, 1.48, 1.42, 1.36, 1.29, 1.23, 1.17, 1.11, 1.04, 0.98},
1029 
1030  {1.82, 1.90, 1.99, 2.07, 2.16, 2.24, 2.33, 2.41, 2.50, 2.58, 2.78, 2.97,
1031  3.16, 3.35, 3.54, 3.74, 3.93, 4.12, 3.75, 3.39, 3.02, 2.66, 2.29, 1.93,
1032  1.56, 1.20, 1.14, 1.09, 1.03, 0.98, 0.92, 0.87, 0.81, 0.76, 0.70},
1033 
1034  {2.10, 2.19, 2.28, 2.38, 2.47, 2.56, 2.65, 2.75, 2.84, 2.93, 3.08, 3.23,
1035  3.38, 3.53, 3.67, 3.82, 3.97, 4.12, 3.71, 3.30, 2.89, 2.49, 2.08, 1.67,
1036  1.26, 0.85, 0.80, 0.75, 0.71, 0.66, 0.61, 0.56, 0.52, 0.47, 0.42},
1037 
1038  {1.82, 1.90, 1.99, 2.07, 2.16, 2.24, 2.33, 2.41, 2.50, 2.58, 2.78, 2.97,
1039  3.16, 3.35, 3.54, 3.74, 3.93, 4.12, 3.75, 3.39, 3.02, 2.66, 2.29, 1.93,
1040  1.56, 1.20, 1.14, 1.09, 1.03, 0.98, 0.92, 0.87, 0.81, 0.76, 0.70},
1041 
1042  {1.54, 1.62, 1.69, 1.77, 1.85, 1.93, 2.00, 2.08, 2.16, 2.24, 2.47, 2.71,
1043  2.94, 3.18, 3.41, 3.65, 3.88, 4.12, 3.80, 3.48, 3.15, 2.83, 2.51, 2.19,
1044  1.87, 1.54, 1.48, 1.42, 1.36, 1.29, 1.23, 1.17, 1.11, 1.04, 0.98},
1045 
1046  {1.26, 1.33, 1.40, 1.47, 1.54, 1.61, 1.68, 1.75, 1.82, 1.89, 2.17, 2.45,
1047  2.73, 3.01, 3.28, 3.56, 3.84, 4.12, 3.84, 3.56, 3.28, 3.01, 2.73, 2.45,
1048  2.17, 1.89, 1.82, 1.75, 1.68, 1.61, 1.54, 1.47, 1.40, 1.33, 1.26},
1049 
1050  {0.98, 1.04, 1.11, 1.17, 1.23, 1.29, 1.36, 1.42, 1.48, 1.54, 1.87, 2.19,
1051  2.51, 2.83, 3.15, 3.48, 3.80, 4.12, 3.88, 3.65, 3.41, 3.18, 2.94, 2.71,
1052  2.47, 2.24, 2.16, 2.08, 2.00, 1.93, 1.85, 1.77, 1.69, 1.62, 1.54},
1053 
1054  {0.70, 0.76, 0.81, 0.87, 0.92, 0.98, 1.03, 1.09, 1.14, 1.20, 1.56, 1.93,
1055  2.29, 2.66, 3.02, 3.39, 3.75, 4.12, 3.93, 3.74, 3.54, 3.35, 3.16, 2.97,
1056  2.78, 2.58, 2.50, 2.41, 2.33, 2.24, 2.16, 2.07, 1.99, 1.90, 1.82}};
1057 
1058  int days[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
1059 
1060  // Set up for the time interpolation:
1061 
1062  if (year % 4 == 0) days[1] = days[1] + 1;
1063 
1064  daymid = days[month] / 2.0;
1065 
1066  if (day >= daymid) {
1067  m1 = month;
1068  m2 = m1 + 1;
1069  if (m2 > 11) m2 = 0;
1070  t = day;
1071  } else {
1072  m2 = month;
1073  m1 = m2 - 1;
1074  if (m1 < 0) m1 = 11;
1075  t = days[m1] + day;
1076  }
1077 
1078  t1 = days[m1] / 2.0;
1079  t2 = days[m1] + (days[m2] / 2.0);
1080 
1081  // Perform the lat. interpolation, and further set up for time interp.:
1082 
1083  if (lat >= 85.0) {
1084  d1 = tabwv[m1][0];
1085  d2 = tabwv[m2][0];
1086  } else if (lat <= (-85.0)) {
1087  d1 = tabwv[m1][34];
1088  d2 = tabwv[m2][34];
1089  } else {
1090  i1 = 16 - (int) (lat / 5.0);
1091  i2 = i1 + 1;
1092  fac = (tabwv[month][i2] - tabwv[month][i1]) / (-5.0);
1093  difflat = lat - (90.0 - (i1 * 5.0));
1094  d1 = tabwv[m1][i1] + (fac * difflat);
1095  d2 = tabwv[m2][i1] + (fac * difflat);
1096  }
1097 
1098  // Complete the time interpolation to get final estimate:
1099 
1100  fac2 = (d2 - d1) / (t2 - t1);
1101 
1102  watvap = d1 + (t - t1) * fac2;
1103 
1104  return watvap;
1105 }
1106 
1107 /********************************************************************************
1108  * Get Earth-Sun distance correction (VarSol)
1109  *
1110  * Calculation of the variability of the solar constant during the
1111  *year. jday is the number of the day in the month dsol is a multiplicative
1112  *factor to apply to the mean value of solar constant.
1113  *
1114  ********************************************************************************/
1115 
1116 float varsol(int32_t jday) {
1117  float om;
1118  float dsol;
1119 
1120  om = (0.9856f * (float) (jday - 4)) * M_PI / 180.f;
1121  dsol = 1.f / powf((1.f - .01673f * cosf(om)), 2.f);
1122 
1123  return dsol;
1124 }
1125 
1126 /*************************************************************************************
1127  * Compute the rise and set time for the given Julian day of year
1128  **
1129  *************************************************************************************/
1130 
1131 void triseset(int32_t jday, float xlon, float xlat, float *trise, float *tset) {
1132  float fac, xlo, xla, xj, a1, a2, et, a3, delta, cosah, ah1, ah2, tsv1, tsv2,
1133  tsm1, tsm2;
1134 
1135  fac = (float) M_PI / 180.f;
1136 
1137  xlo = xlon * fac;
1138  xla = xlat * fac;
1139  xj = (float) (jday - 1);
1140 
1141  a1 = (1.00554f * xj - 6.28306f) * fac;
1142  a2 = (1.93946f * xj + 23.35089f) * fac;
1143  et = -7.67825f * sinf(a1) - 10.09176f * sinf(a2);
1144  a3 = (0.9683f * xj - 78.00878f) * fac;
1145 
1146  delta = 23.4856f * sinf(a3) * fac;
1147  cosah = (-(sinf(xla) * sinf(delta))) / (cosf(xla) * cosf(delta));
1148 
1149  if ((cosah < -1.f) || (cosah > 1.f)) {
1150  *trise = 0.0f;
1151  *tset = 24.0f;
1152  return;
1153  }
1154 
1155  ah1 = acosf(cosah);
1156  ah2 = -ah1;
1157 
1158  tsv1 = ah1 / (15.f * fac);
1159  tsv2 = ah2 / (15.f * fac);
1160  tsm1 = tsv1 + 12.f - et / 60.f;
1161  tsm2 = tsv2 + 12.f - et / 60.f;
1162 
1163  *trise = tsm2 - (xlo / (15.f * fac));
1164  *tset = tsm1 - (xlo / (15.f * fac));
1165  return;
1166 }
1167 
1168 /********************************************************************************
1169  * Get Sun position
1170  ********************************************************************************/
1171 
1172 float get_solz(int32_t jday, float time, float lon, float lat) {
1173  float asol;
1174  double tsm, xla, xj, tet, a1, a2, a3, a4, a5, et, tsv, ah, b1, b2, b3, b4,
1175  b5, b6, b7, delta, amuzero, elev, az, caz, azim, pi2;
1176 
1177  // solar position (zenithal angle asol,azimuthal angle phi0 in degrees)
1178  // j is the day number in the year
1179  //
1180  // mean solar time (heure decimale) (aka decimal time)
1181  double fac = M_PI / 180.;
1182  tsm = time + lon / 15.;
1183  xla = lat * fac;
1184  xj = (float) (jday - 1);
1185  tet = 2. * M_PI * xj / 365.;
1186 
1187  // time equation (in mn.dec)
1188  a1 = 0.000075;
1189  a2 = 0.001868;
1190  a3 = 0.032077;
1191  a4 = 0.014615;
1192  a5 = 0.040849;
1193  et = a1 + a2 * cos(tet) - a3 * sin(tet) - a4 * cos(2. * tet) -
1194  a5 * sin(2. * tet);
1195  et = et * 12. * 60. / M_PI;
1196 
1197  // true solar time
1198 
1199  tsv = tsm + et / 60.;
1200  tsv -= 12.;
1201 
1202  // hour angle
1203 
1204  ah = tsv * 15. * fac;
1205 
1206  // solar declination (in radian)
1207 
1208  b1 = 0.006918;
1209  b2 = 0.399912;
1210  b3 = 0.070257;
1211  b4 = 0.006758;
1212  b5 = 0.000907;
1213  b6 = 0.002697;
1214  b7 = 0 - .001480;
1215  delta = b1 - b2 * cos(tet) + b3 * sin(tet) - b4 * cos(2. * tet) +
1216  b5 * sin(2. * tet) - b6 * cos(3. * tet) + b7 * sin(3. * tet);
1217 
1218  // elevation,azimuth
1219 
1220  amuzero = sin(xla) * sin(delta) + cos(xla) * cos(delta) * cos(ah);
1221  if ((fabs(amuzero) - 1.000) > 0.00000) amuzero = copysign(amuzero, 1.);
1222  elev = asin(amuzero);
1223  az = cos(delta) * sin(ah) / cos(elev);
1224  if ((fabs(az) - 1.000) > 0.00000) az = copysign(az, 1.);
1225 
1226  caz =
1227  (-cos(xla) * sin(delta) + sin(xla) * cos(delta) * cos(ah)) / cos(elev);
1228  azim = asin(az);
1229 
1230  if (caz <= 0.) azim = M_PI - azim;
1231 
1232  if (caz > 0. && az <= 0.) azim = 2. * M_PI + azim;
1233 
1234  azim = azim + M_PI;
1235  pi2 = 2.f * M_PI;
1236 
1237  if (azim > pi2) azim -= pi2;
1238 
1239  elev /= fac;
1240 
1241  // conversion in degrees
1242 
1243  asol = (float) (90. - elev);
1244  // phi0 = azim / fac;
1245 
1246  return asol;
1247 }
1248 
1249 /******************************************************************************
1250  * Given inputs cos(SZ) and tauA500, interpolate the appropriate * surface
1251  *albedo from the table. This is for the case where tau < 1.0 * Inputs are
1252  *assumed to be in the valid ranges (the associated solar * zenith angle
1253  *should be in range 0.0 - 87.134 degrees, and tau * should be in
1254  *the range 0.0 - 0.99)
1255  **
1256  ******************************************************************************/
1257 
1258 float interp_as_taulow(float csz, float tau) {
1259  int s1, s2, t1, t2, i;
1260 
1261  float stot, sdist, ttot, tdist, slope, alb1[2], as;
1262 
1263  float cszt[NS] = {0.05, 0.09, 0.15, 0.21, 0.27, 0.33, 0.39, 0.45,
1264  0.52, 0.60, 0.68, 0.76, 0.84, 0.92, 1.00};
1265 
1266  float taut[NT] = {0.00, 0.05, 0.10, 0.16, 0.24, 0.35, 0.50, 0.70, 0.99};
1267 
1268  // Store the table here rather than read a new file. Fortran stores
1269  // 2-d arrays internally by columns (col1, col2, col3 ...)
1270 
1271  float ast[NS][NT] = {{0.1592253, 0.1258933, 0.1088253, 0.0968594, 0.0880966,
1272  0.0815492, 0.0766266, 0.0721790, 0.0679746},
1273  {0.1944972, 0.1637040, 0.1432702, 0.1240363, 0.1066057,
1274  0.0922638, 0.0826043, 0.0757601, 0.0700167},
1275  {0.1978650, 0.1785707, 0.1643212, 0.1482773, 0.1304572,
1276  0.1119496, 0.0963537, 0.0839010, 0.0740948},
1277  {0.1785221, 0.1673987, 0.1589092, 0.1486409, 0.1359193,
1278  0.1209508, 0.1061981, 0.0920763, 0.0793196},
1279  {0.1531512, 0.1473233, 0.1426633, 0.1366985, 0.1289257,
1280  0.1188203, 0.1078114, 0.0957847, 0.0831087},
1281  {0.1281988, 0.1255739, 0.1234628, 0.1204020, 0.1161201,
1282  0.1101539, 0.1030820, 0.0943791, 0.0841410},
1283  {0.1061354, 0.1052812, 0.1048179, 0.1036617, 0.1016918,
1284  0.0986726, 0.0948040, 0.0894514, 0.0822926},
1285  {0.0877530, 0.0881279, 0.0883850, 0.0884469, 0.0880587,
1286  0.0870923, 0.0854952, 0.0826607, 0.0782380},
1287  {0.0708031, 0.0720173, 0.0727886, 0.0736250, 0.0741367,
1288  0.0746769, 0.0747173, 0.0741294, 0.0722523},
1289  {0.0567974, 0.0582123, 0.0593260, 0.0604172, 0.0615151,
1290  0.0627740, 0.0639047, 0.0647193, 0.0650727},
1291  {0.0472026, 0.0485713, 0.0495830, 0.0507434, 0.0519943,
1292  0.0536504, 0.0551176, 0.0566950, 0.0581553},
1293  {0.0406631, 0.0419177, 0.0429259, 0.0440614, 0.0451823,
1294  0.0467889, 0.0483670, 0.0501810, 0.0522433},
1295  {0.0366926, 0.0377500, 0.0384530, 0.0393431, 0.0404503,
1296  0.0419569, 0.0434430, 0.0451987, 0.0473454},
1297  {0.0343793, 0.0352937, 0.0358116, 0.0364891, 0.0374246,
1298  0.0385732, 0.0398924, 0.0414707, 0.0435983},
1299  {0.0331561, 0.0337733, 0.0341567, 0.0346916, 0.0354239,
1300  0.0364011, 0.0374280, 0.0387133, 0.0405543}};
1301 
1302  // Set up interpolation:
1303 
1304  // Locate subcells in array to use:
1305 
1306  s1 = NS - 2;
1307  s2 = NS - 1;
1308  for (i = 0; i < NS - 1; i++) {
1309  if (csz < cszt[i + 1]) {
1310  s1 = i;
1311  s2 = i + 1;
1312  break;
1313  }
1314  }
1315  t1 = NT - 2;
1316  t2 = NT - 1;
1317  for (i = 0; i < NT - 1; i++) {
1318  if (tau < taut[i + 1]) {
1319  t1 = i;
1320  t2 = i + 1;
1321  break;
1322  }
1323  }
1324 
1325  stot = cszt[s2] - cszt[s1];
1326  sdist = csz - cszt[s1];
1327  ttot = taut[t2] - taut[t1];
1328  tdist = tau - taut[t1];
1329 
1330  slope = (ast[s2][t1] - ast[s1][t1]) / stot;
1331  alb1[0] = ast[s1][t1] + (slope * sdist);
1332  slope = (ast[s2][t2] - ast[s1][t2]) / stot;
1333  alb1[1] = ast[s1][t2] + (slope * sdist);
1334 
1335  slope = (alb1[1] - alb1[0]) / ttot;
1336  as = alb1[0] + (slope * tdist);
1337 
1338  return as;
1339 }
1340 
1341 /********************************************************************************
1342  * Given input cos(SZ), interpolate the appropriate surface albedo
1343  ** from the table. This is for the case where tau > 1.0
1344  ** Input is assumed to be in the valid range (the associated solar
1345  ** zenith angle should be in range 0.0 - 87.134 degrees)
1346  **
1347  ********************************************************************************/
1348 
1349 float interp_as_tauhigh(float csz) {
1350  int s1, s2, i;
1351 
1352  float stot, sdist, slope, as;
1353 
1354  float cszt[NS] = {0.05, 0.09, 0.15, 0.21, 0.27, 0.33, 0.39, 0.45,
1355  0.52, 0.60, 0.68, 0.76, 0.84, 0.92, 1.00};
1356 
1357  float ast[NS] = {0.0623166, 0.0630070, 0.0640739, 0.0650397, 0.0657760,
1358  0.0661690, 0.0659415, 0.0651271, 0.0635101, 0.0610652,
1359  0.0583470, 0.0556146, 0.0530646, 0.0509498, 0.0490149};
1360 
1361  // Set up interpolation:
1362 
1363  // Locate subcells in array to use:
1364 
1365  s1 = NS - 2;
1366  s2 = NS - 1;
1367  for (i = 0; i < NS; i++) {
1368  if (csz < cszt[i + 1]) {
1369  s1 = i;
1370  s2 = i + 1;
1371  break;
1372  }
1373  }
1374 
1375  stot = cszt[s2] - cszt[s1];
1376  sdist = csz - cszt[s1];
1377 
1378  slope = (ast[s2] - ast[s1]) / stot;
1379  as = ast[s1] + (slope * sdist);
1380 
1381  return as;
1382 }
1383 
1384 void getcldalbe(float *TauCld, float *CF, float cosSZ, float t_obs,
1385  float *t_range, float *albe_obs, float *TauCld_obs,
1386  float *CF_obs, size_t t_step, float *wl, size_t bands_num) {
1387  float wave[] = {370., 440., 560., 650., 720.};
1388  float a1[] = {0.78, 0.65, 0.53, 0.49, 0.51};
1389  float b1[] = {0.31, 0.59, 0.79, 0.88, 0.82};
1390  float b2[] = {-0.0782, -0.1212, -0.1054, -0.1276, -0.1054};
1391  float b3[] = {-0.9218, -0.8788, -0.8946, -0.8724, -0.8946};
1392  float k1[] = {1.3, 1.3, 1.3, 1.3, 1.3};
1393  float k2[] = {0.2000, 0.0229, 0.0409, 0.0436, 0.0575};
1394  float k3[] = {0.8000, 1.1754, 1.2436, 1.1836, 1.3637};
1395  float cc[] = {0.0835, 0.1008, 0.1153, 0.1215, 0.1307};
1396  float dd[] = {0.0759, 0.0954, 0.1126, 0.1192, 0.1218};
1397 
1398  size_t s1, s2;
1399 
1400  s1 = t_step - 2;
1401  s2 = t_step - 1;
1402 
1403  float Stot, Sdist, slope;
1404 
1405  search(t_range, 0, t_step - 1, t_obs, &s1, &s2);
1406 
1407  Stot = t_range[s2] - t_range[s1];
1408  Sdist = t_obs - t_range[s1];
1409 
1410  if (t_obs <= t_range[0]) Sdist = 0;
1411  if (t_obs >= t_range[t_step - 1]) Sdist = Stot;
1412  if (s1 != s2)
1413  slope = (TauCld[s2] - TauCld[s1]) / Stot;
1414  else
1415  slope = 0.0;
1416  *TauCld_obs = TauCld[s1] + (slope * Sdist);
1417  if (s1 != s2)
1418  slope = (CF[s2] - CF[s1]) / Stot;
1419  else
1420  slope = 0.0;
1421  *CF_obs = CF[s1] + (slope * Sdist);
1422  // if(isnan(*CF_obs))
1423  // {
1424  // printf("CF_OBS IS NAN,%f, %f, %f ",CF[s1],slope,Sdist);
1425  // exit(EXIT_FAILURE);
1426  // }
1427  float aa, bb;
1428 
1429  float trc[5];
1430 
1431  float alpha = 0; // # alpha - surface albedo
1432 
1433  for (size_t i = 0; i < 5; i++) {
1434  aa = a1[i] + (1 - a1[i]) * exp(-*TauCld_obs * k1[i]);
1435  bb = b1[i] * (1 + b2[i] * exp(-*TauCld_obs * k2[i]) +
1436  b3[i] * exp(-*TauCld_obs * k3[i]));
1437 
1438  trc[i] =
1439  (aa + cosSZ * bb) / (1 + *TauCld_obs * (cc[i] - dd[i] * alpha));
1440  }
1441  for (size_t band = 0; band < bands_num; band++) {
1442  search(wave, 0, 4, wl[band], &s1, &s2);
1443  Stot = wave[s2] - wave[s1];
1444  float Sdist = wl[band] - wave[s1];
1445  float slope = (trc[s2] - trc[s1]) / Stot;
1446  if (s1 == s2) slope = 0.0;
1447  float trc_obs = trc[s1] + (slope * Sdist);
1448  albe_obs[band] = (1 - trc_obs) * (*CF_obs);
1449  }
1450 }
1451 
1452 // ##############################################################################
1453 // # new parameterization of As, from Jin et al. (2011)
1454 // # Jin, Z., Y. Qiao, Y. Wang, Y. Fang, and W. Yi, 2011: A new parameterization
1455 // of # spectral and broadband ocean surface albedo. Opt. Express, 19, no. 27,
1456 // 26429-26443, # doi:10.1364/OE.19.026429.
1457 // ##############################################################################
1458 
1459 // # Regression function (Equation 4)
1460 
1461 float ff(float u, float v) {
1462  float p[] = {0.0152, -1.7873, 6.8972, -8.5778, 4.071, -7.6446,
1463  0.1643, -7.8409, -3.5639, -2.3588, 10.0538};
1464  const float zmod =
1465  (p[0] + p[1] * u + p[2] * pow(u, 2) + p[3] * pow(u, 3) + p[4] * v +
1466  p[5] * u * v) *
1467  exp(p[6] + p[7] * u + p[8] * pow(u, 2) + p[9] * v + p[10] * u * v);
1468  return zmod;
1469 }
1470 
1475 float getrefm(float wl, const luts_par *luts_data) {
1476  size_t dims[] = {luts_data->soa_lut.dim_wlt_reft};
1477  float *grid[1];
1478  grid[0] = luts_data->soa_lut.wlt_reft;
1479  float point[] = {wl};
1480  float refm = interp1d(dims, point, grid, luts_data->soa_lut.reft);
1481  if (refm > 1.34) refm = 1.34;
1482  return refm;
1483 }
1484 
1485 float getr0(float wl, float chl, float u0, const luts_par *luts_data) {
1486  float *grid[1];
1487  float point[] = {wl};
1488  float chlabs, bw, aw;
1489  // get chlor
1490  {
1491  size_t dims[] = {luts_data->soa_lut.dim_jwvl};
1492  grid[0] = luts_data->soa_lut.jwl;
1493  chlabs = interp1d(dims, point, grid, luts_data->soa_lut.achl);
1494  }
1495 
1496  // get bw
1497  {
1498  size_t dims[] = {luts_data->soa_lut.dim_wv1};
1499  grid[0] = luts_data->soa_lut.wv1;
1500  bw = interp1d(dims, point, grid, luts_data->soa_lut.bw1);
1501  }
1502 
1503  // get aw
1504  {
1505  size_t dims[] = {luts_data->soa_lut.dim_wv3};
1506  grid[0] = luts_data->soa_lut.wv3;
1507  aw = interp1d(dims, point, grid, luts_data->soa_lut.aw3);
1508  }
1509  const float ylmd = exp(0.014 * (440.0 - wl));
1510  const float aw440 = 0.00635;
1511  const float ap = 0.06 * chlabs * pow(chl, 0.65) +
1512  0.2 * (aw440 + 0.06 * pow(chl, 0.65)) * ylmd;
1513 
1514  const float bp550 = 0.416 * pow(chl, 0.766);
1515 
1516  float v, bbp;
1517 
1518  if (chl >= 2) {
1519  v = 0.;
1520  bbp = (0.002 + 0.01 * (0.5 - 0.25 * log10(chl))) * bp550;
1521  } else {
1522  if (chl > 0.02) {
1523  v = 0.5 * (log10(chl) - 0.3);
1524  bbp =
1525  (0.002 + 0.01 * (0.5 - 0.25 * log10(chl)) * pow(wl / 550., v)) *
1526  bp550;
1527  } else {
1528  bbp = 0.019 * (550. / wl) * bp550;
1529  }
1530  }
1531 
1532  // # hb=h/(h+2*bbpf*(1.-h)) ;Morel-Gentili(1991), Eq (12)
1533  float hb = 0.5 * bw / (0.5 * bw + bbp);
1534 
1535  // #; use Morel 91 to get f
1536 
1537  float f =
1538  0.6279 - 0.2227 * hb - 0.0513 * hb * hb + (-0.3119 + 0.2465 * hb) * u0;
1539 
1540  float r0 = f * (0.5 * bw + bbp) / (aw + ap);
1541  return r0;
1542 }
1543 
1544 float getosa(float wl, float sza, float wind, float chl, float fr,
1545  const luts_par *luts_data) { // # wl wavelength (nm), set wl=0 for broadband
1546  // # chl chlorophyll (g/m3)
1547  // # sza solar zenith angle (degree)
1548  // # wind wind speed (m/s)
1549  // # fr fraction of diffuse incidence
1550  const float sig = sqrt(0.003 + 0.00512 * wind); // #convert to Sigma
1551 
1552  const float refm = getrefm(wl, luts_data);
1553 
1554  const float csz = cos(sza * M_PI / 180.);
1555 
1556  const float xx2 = sqrt(1.0 - (1.0 - csz * csz) / refm / refm);
1557 
1558  float rr0 = (0.5 * (pow(((xx2 - csz * refm) / (xx2 + csz * refm)), 2) +
1559  pow(((csz - refm * xx2) / (csz + refm * xx2)), 2)));
1560 
1561  float rrr = 0.5 * (pow(((xx2 - 1.34 * csz) / (xx2 + 1.34 * csz)), 2) +
1562  pow(((csz - 1.34 * xx2) / (csz + 1.34 * xx2)), 2));
1563 
1564  const float r11 = rr0 - ff(csz, sig) * rr0 / rrr; // # direct surface
1565  // albedo
1566 
1567  // # Water volume scattering
1568  const float r22 = 0.48168549 - 0.014894708 * sig - 0.20703885 * sig * sig;
1569 
1570  float r00 = getr0(wl, chl, csz, luts_data);
1571 
1572  const float rw = r00 * (1. - r22) * (1. - r11) /
1573  (1. - r00 * r22); // #direct water-leaving albedo
1574 
1575  const float ue =
1576  0.676; // # the equivalent u_unif for diffuse incidence
1577  const float ue2 = sqrt(1.0 - (1.0 - ue * ue) / refm / refm);
1578 
1579  rr0 = 0.5 * (pow(((ue2 - refm * ue) / (ue2 + refm * ue)), 2) +
1580  pow(((ue - refm * ue2) / (ue + refm * ue2)), 2));
1581  rrr = 0.5 * (pow(((ue2 - 1.34 * ue) / (ue2 + 1.34 * ue)), 2) +
1582  pow(((ue - 1.34 * ue2) / (ue + 1.34 * ue2)), 2));
1583 
1584  float r11df = rr0 - ff(ue, sig) * rr0 / rrr;
1585 
1586  r00 = getr0(wl, chl, ue, luts_data);
1587 
1588  float rwdf = r00 * (1. - r22) * (1. - r11df) /
1589  (1. - r00 * r22); // # diffuse water-leaving albedo
1590 
1591  float pc[] = {-0.1482, -0.012, 0.1609, -0.0244};
1592 
1593  float rdf;
1594  if (fr > 0.99)
1595  rdf = -0.1479 + 0.1502 * refm - 0.0176 * sig * refm;
1596  else
1597  rdf = pc[0] + pc[1] * sig + pc[2] * refm +
1598  pc[3] * sig * refm; // # surface diffuse (Eq 5a-5b)
1599 
1600  // float rwt =
1601  // (1. - fr) * rw +
1602  // fr * rwdf; // #total parameterized water-leaving albedo
1603 
1604  float albt = (1. - fr) * (r11 + rw) +
1605  fr * (rdf + rwdf); // # total parameterized albedo
1606 
1607  // # Foam correction (Eq 16-17)
1608  float fwc = (2.95e-6) * pow(wind, 3.52);
1609 
1610  float osa = (1. - fwc) * albt + fwc * 0.55;
1611  return osa;
1612 }
1613 
1614 float SunGlint(float sz, float vz, float ra, float ws) {
1615  const float ssz = sin(sz * M_PI / 180.);
1616  const float svz = sin(vz * M_PI / 180.);
1617  const float csz = cos(sz * M_PI / 180.);
1618  const float cvz = cos(vz * M_PI / 180.);
1619  const float cra = cos(ra * M_PI / 180.);
1620 
1621  // # Determine omega
1622  const float cos_2omega = csz * cvz + ssz * svz * cra;
1623  const float omega = acos(cos_2omega) / 2.0; // # in radians
1624  const float sOmega = sin(omega);
1625  const float cOmega = cos(omega);
1626  const float omegaP = asin(sOmega / 1.34); // # in radians, Snell's law
1627  const float omegaDp = omega + omegaP;
1628  const float omegaDm = omega - omegaP;
1629 
1630  const float term1 = sin(omegaDm) / sin(omegaDp);
1631  const float term2 = tan(omegaDm) / tan(omegaDp);
1632  float r_omega = 0.5 * (term1 * term1 + term2 * term2);
1633  if (omega < 0.0174533) r_omega = 0.0211119;
1634 
1635  const float beta = acos((csz + cvz) / (2.0 * cOmega)); // # in radians
1636  const float cBeta = cos(beta);
1637  const float tBeta = tan(beta);
1638  const float sig2 = 0.003 + 0.00512 * ws;
1639  const float P_V_beta =
1640  (1.0 / (2.0 * sig2 * M_PI)) * exp(-tBeta * tBeta / (2.0 * sig2));
1641  const float rho_g =
1642  M_PI * P_V_beta * r_omega / (4.0 * csz * cvz * pow((cBeta), 4));
1643  return rho_g;
1644 }
1645 
1646 void set_searh_grid(size_t pts[][2], float width[], size_t ndims,
1647  const size_t *dims, const float *point, float **grid) {
1648  for (size_t i_dim = 0; i_dim < ndims; i_dim++) {
1649  size_t dim_size = dims[i_dim];
1650  size_t st = 0;
1651  size_t end = dim_size - 1;
1652  search(grid[i_dim], st, end, point[i_dim], &st, &end);
1653  pts[i_dim][0] = st;
1654  pts[i_dim][1] = end;
1655  }
1656  for (size_t i_dim = 0; i_dim < ndims; i_dim++) {
1657  size_t st = pts[i_dim][0];
1658  size_t end = pts[i_dim][1];
1659  if (st != end)
1660  width[i_dim] = (point[i_dim] - grid[i_dim][st]) /
1661  (grid[i_dim][end] - grid[i_dim][st]);
1662  else
1663  width[i_dim] = 0.0;
1664  // { std::cout << "width[i_dim]" << width[i_dim] << std::endl; }
1665  }
1666 }
1667 
1668 size_t get6dindex(const size_t *point, const size_t *dims) {
1669  return point[0] * dims[1] * dims[2] * dims[3] * dims[4] * dims[5] +
1670  point[1] * dims[2] * dims[3] * dims[4] * dims[5] +
1671  point[2] * dims[3] * dims[4] * dims[5] +
1672  point[3] * dims[4] * dims[5] + point[4] * dims[5] + point[5];
1673 }
1674 
1675 float interp6d(const size_t *dims, const float *point, float **grid,
1676  const float *lut) {
1677  const size_t ndims = 6;
1678  size_t pts[ndims][2];
1679  float width[ndims];
1680  set_searh_grid(pts, width, ndims, dims, point, grid);
1681  float temp0[2][2][2][2][2][2];
1682  for (size_t i0 = 0; i0 < 2; i0++)
1683  for (size_t i1 = 0; i1 < 2; i1++)
1684  for (size_t i2 = 0; i2 < 2; i2++)
1685  for (size_t i3 = 0; i3 < 2; i3++)
1686  for (size_t i4 = 0; i4 < 2; i4++)
1687  for (size_t i5 = 0; i5 < 2; i5++) {
1688  size_t point_pos[] = {pts[0][i0], pts[1][i1],
1689  pts[2][i2], pts[3][i3],
1690  pts[4][i4], pts[5][i5]};
1691  temp0[i0][i1][i2][i3][i4][i5] =
1692  lut[get6dindex(point_pos, dims)];
1693  }
1694  float temp1[2][2][2][2][2];
1695  for (size_t i1 = 0; i1 < 2; i1++)
1696  for (size_t i2 = 0; i2 < 2; i2++)
1697  for (size_t i3 = 0; i3 < 2; i3++)
1698  for (size_t i4 = 0; i4 < 2; i4++)
1699  for (size_t i5 = 0; i5 < 2; i5++) {
1700  temp1[i1][i2][i3][i4][i5] =
1701  (temp0[1][i1][i2][i3][i4][i5] -
1702  temp0[0][i1][i2][i3][i4][i5]) *
1703  width[0] +
1704  temp0[0][i1][i2][i3][i4][i5];
1705  }
1706 
1707  float temp2[2][2][2][2];
1708  for (size_t i2 = 0; i2 < 2; i2++)
1709  for (size_t i3 = 0; i3 < 2; i3++)
1710  for (size_t i4 = 0; i4 < 2; i4++)
1711  for (size_t i5 = 0; i5 < 2; i5++) {
1712  temp2[i2][i3][i4][i5] =
1713  (temp1[1][i2][i3][i4][i5] - temp1[0][i2][i3][i4][i5]) *
1714  width[1] +
1715  temp1[0][i2][i3][i4][i5];
1716  }
1717 
1718  float temp3[2][2][2];
1719  for (size_t i3 = 0; i3 < 2; i3++)
1720  for (size_t i4 = 0; i4 < 2; i4++)
1721  for (size_t i5 = 0; i5 < 2; i5++) {
1722  temp3[i3][i4][i5] =
1723  (temp2[1][i3][i4][i5] - temp2[0][i3][i4][i5]) * width[2] +
1724  temp2[0][i3][i4][i5];
1725  }
1726  float temp4[2][2];
1727  for (size_t i4 = 0; i4 < 2; i4++)
1728  for (size_t i5 = 0; i5 < 2; i5++) {
1729  temp4[i4][i5] = (temp3[1][i4][i5] - temp3[0][i4][i5]) * width[3] +
1730  temp3[0][i4][i5];
1731  }
1732 
1733  float temp5[2];
1734  for (size_t i5 = 0; i5 < 2; i5++) {
1735  temp5[i5] = (temp4[1][i5] - temp4[0][i5]) * width[4] + temp4[0][i5];
1736  }
1737 
1738  float temp6;
1739  temp6 = (temp5[1] - temp5[0]) * width[5] + temp5[0];
1740 
1741  return temp6;
1742 }
1743 
1744 size_t get4dindex(const size_t *point, const size_t *dims) {
1745  return point[0] * dims[1] * dims[2] * dims[3] +
1746  point[1] * dims[2] * dims[3] + point[2] * dims[3] + point[3];
1747 }
1748 
1749 float interp4d(const size_t *dims, const float *point, float **grid,
1750  const float *lut) {
1751  const size_t ndims = 4;
1752  size_t pts[ndims][2];
1753  float width[ndims];
1754  set_searh_grid(pts, width, ndims, dims, point, grid);
1755  float temp0[2][2][2][2];
1756  for (size_t i0 = 0; i0 < 2; i0++)
1757  for (size_t i1 = 0; i1 < 2; i1++)
1758  for (size_t i2 = 0; i2 < 2; i2++)
1759  for (size_t i3 = 0; i3 < 2; i3++) {
1760  size_t point_pos[] = {pts[0][i0], pts[1][i1], pts[2][i2],
1761  pts[3][i3]};
1762  temp0[i0][i1][i2][i3] = lut[get4dindex(point_pos, dims)];
1763  }
1764  float temp1[2][2][2];
1765  for (size_t i1 = 0; i1 < 2; i1++)
1766  for (size_t i2 = 0; i2 < 2; i2++)
1767  for (size_t i3 = 0; i3 < 2; i3++) {
1768  temp1[i1][i2][i3] =
1769  (temp0[1][i1][i2][i3] - temp0[0][i1][i2][i3]) * width[0] +
1770  temp0[0][i1][i2][i3];
1771  }
1772 
1773  float temp2[2][2];
1774  for (size_t i2 = 0; i2 < 2; i2++)
1775  for (size_t i3 = 0; i3 < 2; i3++) {
1776  temp2[i2][i3] = (temp1[1][i2][i3] - temp1[0][i2][i3]) * width[1] +
1777  temp1[0][i2][i3];
1778  }
1779 
1780  float temp3[2];
1781  for (size_t i3 = 0; i3 < 2; i3++) {
1782  temp3[i3] = (temp2[1][i3] - temp2[0][i3]) * width[2] + temp2[0][i3];
1783  }
1784 
1785  float temp4;
1786  temp4 = (temp3[1] - temp3[0]) * width[3] + temp3[0];
1787 
1788  return temp4;
1789 }
1790 
1791 size_t get3dindex(const size_t *point, const size_t *dims) {
1792  return point[0] * dims[1] * dims[2] + point[1] * dims[2] + point[2];
1793 }
1794 
1795 float interp3d(const size_t *dims, const float *point, float **grid,
1796  const float *lut) {
1797  const size_t ndims = 3;
1798  size_t pts[ndims][2];
1799  float width[ndims];
1800  set_searh_grid(pts, width, ndims, dims, point, grid);
1801  float temp0[2][2][2];
1802  for (size_t i0 = 0; i0 < 2; i0++)
1803  for (size_t i1 = 0; i1 < 2; i1++)
1804  for (size_t i2 = 0; i2 < 2; i2++) {
1805  size_t point_pos[] = {pts[0][i0], pts[1][i1], pts[2][i2]};
1806  temp0[i0][i1][i2] = lut[get3dindex(point_pos, dims)];
1807  }
1808  float temp1[2][2];
1809  for (size_t i1 = 0; i1 < 2; i1++)
1810  for (size_t i2 = 0; i2 < 2; i2++) {
1811  temp1[i1][i2] = (temp0[1][i1][i2] - temp0[0][i1][i2]) * width[0] +
1812  temp0[0][i1][i2];
1813  }
1814 
1815  float temp2[2];
1816  for (size_t i2 = 0; i2 < 2; i2++) {
1817  temp2[i2] = (temp1[1][i2] - temp1[0][i2]) * width[1] + temp1[0][i2];
1818  }
1819 
1820  float temp3 = (temp2[1] - temp2[0]) * width[2] + temp2[0];
1821 
1822  return temp3;
1823 }
1824 
1825 size_t get1dindex(const size_t *point, const size_t *dims) { return point[0]; }
1826 
1827 float interp1d(const size_t *dims, const float *point, float **grid,
1828  const float *lut) {
1829  const size_t ndims = 1;
1830  size_t pts[ndims][2];
1831  float width[ndims];
1832  set_searh_grid(pts, width, ndims, dims, point, grid);
1833  float temp0[2];
1834  for (size_t i0 = 0; i0 < 2; i0++) {
1835  size_t point_pos[] = {pts[0][i0]};
1836  temp0[i0] = lut[get1dindex(point_pos, dims)];
1837  }
1838 
1839  float temp1;
1840  temp1 = (temp0[1] - temp0[0]) * width[0] + temp0[0];
1841 
1842  return temp1;
1843 }
hyper_cube_linear_interp
float hyper_cube_linear_interp(size_t n_dims, float *hyper_cube, float *width)
Definition: par_utils.c:131
sensorId2SensorDir
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:315
t
data_t t[NROOTS+1]
Definition: decode_rs.h:77
get4dindex
size_t get4dindex(const size_t *point, const size_t *dims)
Definition: par_utils.c:1744
get1dindex
size_t get1dindex(const size_t *point, const size_t *dims)
Definition: par_utils.c:1825
j
int j
Definition: decode_rs.h:73
day
int32_t day
Definition: atrem_corl1v2.h:308
status
int status
Definition: l1_czcs_hdf.c:32
NPHASE
#define NPHASE
Definition: dtranbrdf.cpp:10
NT
#define NT
Definition: par_utils.c:45
NMODEL
#define NMODEL
Definition: par_utils.c:43
openDS
idDS openDS(const char *filename)
Definition: wrapper.c:606
allocateMemory
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
FALSE
#define FALSE
Definition: rice.h:164
FILE_ERROR
#define FILE_ERROR(str, name)
Definition: par_utils.c:53
NULL
#define NULL
Definition: decode_rs.h:63
interp3d
float interp3d(const size_t *dims, const float *point, float **grid, const float *lut)
3-dimensional interpolation
Definition: par_utils.c:1795
interp6d
float interp6d(const size_t *dims, const float *point, float **grid, const float *lut)
6-dimensional interpolation
Definition: par_utils.c:1675
band
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second band
Definition: HOWTO_Add_a_product.txt:42
color_dtdb.v
v
Definition: color_dtdb.py:430
core_arrays::latitude
real(single), dimension(:,:), allocatable latitude
Definition: core_arrays.f90:149
unix2ymds
void unix2ymds(double usec, int16_t *year, int16_t *mon, int16_t *day, double *secs)
Definition: unix2ymds.c:8
kasten_equation
float kasten_equation(float solz)
Kasten eq to compute airmass.
Definition: par_utils.c:58
get_hyper_cube_index
size_t get_hyper_cube_index(size_t n_dims, const size_t *dims, const size_t *indexes)
Definition: par_utils.c:121
read_aerosol_par
void read_aerosol_par(l2str *l2rec, int32_t nbands, float *tablewavelengths, float *tablephaseangles, float *tablealphas, float *tableomegas, float *tableaerphasefunc)
Definition: par_utils.c:802
wavelength
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending wavelength
Definition: HOWTO_Add_a_product.txt:24
TRUE
#define TRUE
Definition: rice.h:165
as
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
jday
int32_t jday(int16_t i, int16_t j, int16_t k)
Converts a calendar date to the corresponding Julian day starting at noon on the calendar date....
Definition: jday.c:14
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
get6dindex
size_t get6dindex(const size_t *point, const size_t *dims)
Definition: par_utils.c:1668
GetAerPhase
void GetAerPhase(l2str *l2rec, int ip, int32_t nbands, float angstrom, float *aerphasefunc, float *omega, float *modelAngstrom)
Definition: par_utils.c:680
M_PI
#define M_PI
Definition: dtranbrdf.cpp:19
fac
#define fac
Definition: atmocor1_land.c:116
getrefm
float getrefm(float wl, const luts_par *luts_data)
Definition: par_utils.c:1475
ff
float ff(float u, float v)
Definition: par_utils.c:1461
run_local.lat
lat
Definition: run_local.py:69
f
double precision function f(R1)
Definition: tmd.lp.f:1454
core_arrays::solar_zenith_angle
real, dimension(:,:), allocatable solar_zenith_angle
Definition: core_arrays.f90:6
getosa
float getosa(float wl, float sza, float wind, float chl, float fr, const luts_par *luts_data)
Definition: par_utils.c:1544
get_luts_data
void get_luts_data(l2str *l2rec, luts_par *luts_data)
Get the luts data object For now it includes the MERRA data. Should be moved after the merra data is ...
Definition: par_utils.c:244
run_local.lon
lon
Definition: run_local.py:69
interp_as_tauhigh
float interp_as_tauhigh(float csz)
Definition: par_utils.c:1349
get_nc_var
void get_nc_var(const idDS *struct_nc, const char *name, float *var, const size_t *start, const size_t *count)
Definition: par_utils.c:227
omega
data_t omega[NROOTS+1]
Definition: decode_rs.h:77
getcldalbe
void getcldalbe(float *TauCld, float *CF, float cosSZ, float t_obs, float *t_range, float *albe_obs, float *TauCld_obs, float *CF_obs, size_t t_step, float *wl, size_t bands_num)
//
Definition: par_utils.c:1384
search
size_t search(const float *arr, size_t s, size_t e, float val, size_t *i_s, size_t *i_e)
Binary search algorithm.
Definition: par_utils.c:64
slope
float32 slope[]
Definition: l2lists.h:30
delta
const double delta
Definition: VcstCmnConsts.h:523
get_nc_dim
void get_nc_dim(const idDS *struct_nc, const char *name, size_t *len)
Definition: par_utils.c:217
unix2yds
void unix2yds(double usec, short *year, short *day, double *secs)
luts_data
luts_par luts_data
Definition: calc_par.c:89
interp1d
float interp1d(const size_t *dims, const float *point, float **grid, const float *lut)
1-dimensional interpolation
Definition: par_utils.c:1827
color_dtdb.result
result
Definition: color_dtdb.py:197
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
color_dtdb.st
list st
Definition: color_dtdb.py:273
b1
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
NS
#define NS
Definition: par_utils.c:44
getr0
float getr0(float wl, float chl, float u0, const luts_par *luts_data)
Definition: par_utils.c:1485
xlon
float xlon[LAC_PIXEL_NUM]
Definition: l1a_seawifs.c:93
triseset
void triseset(int32_t jday, float xlon, float xlat, float *trise, float *tset)
Definition: par_utils.c:1131
get_reduced_product
size_t get_reduced_product(size_t n_dims, const size_t *dims)
Definition: par_utils.c:113
EstimateWatVap
float EstimateWatVap(int32_t year, int32_t month, int32_t day, float lat)
Definition: par_utils.c:1002
core_arrays::relative_azimuth_angle
real(single), dimension(:,:), allocatable relative_azimuth_angle
Definition: core_arrays.f90:118
interp4d
float interp4d(const size_t *dims, const float *point, float **grid, const float *lut)
4-dimensional interpolation
Definition: par_utils.c:1749
get_solz
float get_solz(int32_t jday, float time, float lon, float lat)
Definition: par_utils.c:1172
beta
double beta
Definition: GEO_earth_location.c:14
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
u
data_t u
Definition: decode_rs.h:74
time
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
fabs
#define fabs(a)
Definition: misc.h:93
get3dindex
size_t get3dindex(const size_t *point, const size_t *dims)
Definition: par_utils.c:1791
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
idDS
Definition: dfutils.h:28
EstimateDobson
float EstimateDobson(int32_t year, int32_t month, int32_t day, float lat)
Definition: par_utils.c:876
s
data_t s[NROOTS]
Definition: decode_rs.h:75
set_searh_grid
void set_searh_grid(size_t pts[][2], float width[], size_t ndims, const size_t *dims, const float *point, float **grid)
Definition: par_utils.c:1646
color_dtdb.alpha
alpha
Definition: color_dtdb.py:352
interpnd
float interpnd(size_t n_dims, const size_t *dims, const float *point, float **grid, const float *lut)
N-dimensional interpolation.
Definition: par_utils.c:158
check_products.end
end
Definition: check_products.py:189
endDS
int endDS(idDS ds_id)
Definition: wrapper.c:624
i
int i
Definition: decode_rs.h:71
viirs_ler_luts::p1
real p1
Definition: viirs_ler_luts.f95:37
SunGlint
float SunGlint(float sz, float vz, float ra, float ws)
Definition: par_utils.c:1614
varsol
float varsol(int32_t jday)
Definition: par_utils.c:1116
par_utils.h
viirs_ler_luts::sza
real sza
Definition: viirs_ler_luts.f95:46
interp_as_taulow
float interp_as_taulow(float csz, float tau)
Definition: par_utils.c:1258
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
viirs_ler_luts::xlat
real xlat
Definition: viirs_ler_luts.f95:44
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

Mission Quality Monitor

HELP

How to Use Earthdata Search

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

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster