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
get_Kd.c
Go to the documentation of this file.
1 /* =================================================================== */
2 /* Algorithms for computing diffuse attenuation coefficient for MSl12. */
3 /* */
4 /* B. Franz, NASA Ocean Biology Processing Group, SAIC, March 2005. */
5 /* =================================================================== */
6 
7 #include <stdlib.h>
8 #include <math.h>
9 #include "l12_proto.h"
10 
11 #define KD_MAX 6.4
12 #define KD_MIN 0.016
13 
14 static float kdbad = BAD_FLT;
15 static float *Kd_spectral;
16 
17 /* ------------------------------------------------------------------- */
18 /* Kd490_KD2 - diffuse attenuation at 490nm (2-band polynomial). */
19 /* */
20 /* Inputs: */
21 /* l2rec - level-2 structure containing one complete scan after */
22 /* atmospheric correction. */
23 /* Outputs: */
24 /* k490 - diffuse attenuation coefficient, 1 value per pixel. */
25 /* */
26 /* Algorithm: P.J.Werdell, June 2009 */
27 
28 /* ------------------------------------------------------------------- */
29 void Kd490_KD2(l2str *l2rec, float *Kd) {
30  static int32_t *w = NULL;
31  static float *a = NULL;
32  static int ib1 = -1;
33  static int ib2 = -1;
34 
35  int32_t ip, ipb;
36  float R;
37 
38  l1str *l1rec = l2rec->l1rec;
39 
40  if (w == NULL) {
41  w = input->kd2w;
42  a = input->kd2c;
43  if (w[0] < 0 || w[1] < 0) {
44  printf("Kd490_KD2: algorithm coefficients not provided for this sensor.\n");
45  exit(1);
46  }
47  ib1 = bindex_get(w[0]);
48  ib2 = bindex_get(w[1]);
49  if (ib1 < 0 || ib2 < 0) {
50  printf("Kd490_KD2: incompatible sensor wavelengths for this algorithm\n");
51  exit(1);
52  }
53  }
54 
55  for (ip = 0; ip < l1rec->npix; ip++) {
56 
57  ipb = l1rec->l1file->nbands*ip;
58 
59  if (l1rec->mask[ip] ||
60  l2rec->Rrs[ipb + ib1] <= 0.0 || l2rec->Rrs[ipb + ib2] <= 0.0) {
61  Kd[ip] = kdbad;
62  l1rec->flags[ip] |= PRODFAIL;
63  } else {
64  R = log10(l2rec->Rrs[ipb + ib1] / l2rec->Rrs[ipb + ib2]);
65  if (isnan(R)) {
66  Kd[ip] = kdbad;
67  l1rec->flags[ip] |= PRODFAIL;
68  } else {
69  Kd[ip] = a[0] + pow(10.0, a[1] + R * (a[2] + R * (a[3] + R * (a[4] + R * a[5]))));
70  if (Kd[ip] > KD_MAX) {
71  Kd[ip] = KD_MAX;
72  l1rec->flags[ip] |= PRODWARN;
73  }
74  }
75  }
76  }
77 }
78 
79 /* ------------------------------------------------------------------- */
80 /* unc_Kd490_KD2 - uncertainty of diffuse attenuation at 490nm */
81 /* (2-band polynomial). */
82 /* */
83 /* Inputs: */
84 /* l2rec - level-2 structure containing one complete scan after */
85 /* atmospheric correction. */
86 /* Outputs: */
87 /* uk490 - standard uncertainty in the diffuse attenuation */
88 /* coefficientat 490 nm, 1 value per pixel. Model uses */
89 /* first order analytical propagation */
90 /* */
91 /* Algorithm: P.J.Werdell, June 2009 */
92 /* Uncertinay method: McKinna et al (2019) */
93 
94 /* ------------------------------------------------------------------- */
95 void unc_Kd490_KD2(l2str *l2rec, float *uKd) {
96  static int32_t *w = NULL;
97  static float *a = NULL;
98  static int ib1 = -1;
99  static int ib2 = -1;
100 
101  int32_t ip, ipb;
102  float R, Kd, X;
103  float dRdR1, dRdR2;
104  float dXdR, dKdX;
105  float dKdR1, dKdR2;
106  float Rrs1, Rrs2;
107  float uRrs1 =0.0;
108  float uRrs2 = 0.0;
109 
110  l1str *l1rec = l2rec->l1rec;
111  uncertainty_t *uncertainty=l1rec->uncertainty;
112 
113  if (w == NULL) {
114  w = input->kd2w;
115  a = input->kd2c;
116  if (w[0] < 0 || w[1] < 0) {
117  printf("Kd490_KD2: algorithm coefficients not provided for this sensor.\n");
118  exit(1);
119  }
120  ib1 = bindex_get(w[0]);
121  ib2 = bindex_get(w[1]);
122  if (ib1 < 0 || ib2 < 0) {
123  printf("Kd490_KD2: incompatible sensor wavelengths for this algorithm\n");
124  exit(1);
125  }
126  }
127 
128 
129  for (ip = 0; ip < l1rec->npix; ip++) {
130 
131  ipb = l1rec->l1file->nbands*ip;
132 
133  Rrs1 = l2rec->Rrs[ipb + ib1];
134  Rrs2 = l2rec->Rrs[ipb + ib2];
135 
136  if (uncertainty) {
137  uRrs1 = l2rec->Rrs[ipb + ib1];
138  uRrs2 = l2rec->Rrs[ipb + ib2];
139  }
140 
141  if (l1rec->mask[ip] || Rrs1 <= 0.0 || Rrs2 <= 0.0) {
142  uKd[ip] = kdbad;
143  l1rec->flags[ip] |= PRODFAIL;
144  } else {
145 
146  R = log10(Rrs1 / Rrs2);
147  dRdR1 = 1./(log(10)*Rrs1);
148  dRdR2 = -1./(log(10)*Rrs2);
149 
150  if (isnan(R)) {
151  uKd[ip] = kdbad;
152  l1rec->flags[ip] |= PRODFAIL;
153  } else {
154  X = a[1] + R * (a[2] + R * (a[3] + R * (a[4] + R * a[5])));
155  dXdR = a[2] + R * (2*a[3] + R*(3*a[4] + R*(4*a[5])));
156  Kd = a[0] + pow(10.0, X);
157  dKdX = log(10)*pow(10.,X);
158  dKdR1 = dKdX * dXdR * dRdR1;
159  dKdR2 = dKdX * dXdR * dRdR2;
160  uKd[ip] = sqrt( pow(dKdR1* uRrs1,2.) + pow(dKdR2*uRrs2,2.) );
161 
162  if (Kd > KD_MAX) {
163  uKd[ip] = 0.0;
164  l1rec->flags[ip] |= PRODWARN;
165  }
166  }
167  }
168  }
169 }
170 
171 /* ------------------------------------------------------------------- */
172 /* Kd490_mueller() - diffuse attenuation at 490nm (J. Mueller). */
173 /* */
174 /* Inputs: */
175 /* l2rec - level-2 structure containing one complete scan after */
176 /* atmospheric correction. */
177 /* Outputs: */
178 /* k490 - diffuse attenuation coefficient, 1 value per pixel. */
179 /* */
180 /* Algorithm provided by: J. Mueller */
181 /* OCTS coefficents: S. Bailey, 16 July 2001 */
182 
183 /* ------------------------------------------------------------------- */
184 void Kd490_mueller(l2str *l2rec, float k490[]) {
185  int32_t ip, ipb;
186 
187  static float a = 0.15645;
188  static float b = -1.5401;
189  static float Kw = 0.016;
190 
191  static float maxval = 6.4;
192 
193  l1str *l1rec = l2rec->l1rec;
194 
195  if (l1rec->l1file->sensorID == OCTS) {
196  /* Coefs for 490/565 band combination */
197  a = 0.2166;
198  b = -1.6355;
199  }
200 
201  for (ip = 0; ip < l1rec->npix; ip++) {
202 
203  ipb = l1rec->l1file->nbands*ip;
204 
205  if (l1rec->mask[ip]) { /* pixel was masked */
206  k490[ip] = kdbad;
207  l1rec->flags[ip] |= PRODFAIL;
208  } else if (l2rec->nLw[ipb + 2] <= 0.0) { /* unknown, high attenuation */
209  k490[ip] = maxval;
210  } else if (l2rec->nLw[ipb + 4] <= 0.0) { /* unknown, low attenuation */
211  k490[ip] = Kw;
212  } else {
213  k490[ip] = a * pow(l2rec->nLw[ipb + 2] / l2rec->nLw[ipb + 4], b) + Kw;
214  if (k490[ip] > maxval) {
215  k490[ip] = maxval;
216  }
217  }
218  }
219 }
220 
221 
222 /* ------------------------------------------------------------------- */
223 /* Kd490_obpg - diffuse attenuation at 490nm (Mueller & Werdell). */
224 /* */
225 /* Inputs: */
226 /* l2rec - level-2 structure containing one complete scan after */
227 /* atmospheric correction. */
228 /* Outputs: */
229 /* k490 - diffuse attenuation coefficient, 1 value per pixel. */
230 /* */
231 /* Algorithm provided by: J. Mueller */
232 /* New coefficents and updated form: P.J.Werdell, February 2005. */
233 
234 /* ------------------------------------------------------------------- */
235 void Kd490_obpg(l2str *l2rec, float k490[]) {
236  int32_t ip, ipb;
237 
238  static float a;
239  static float b;
240  static int32_t ib490;
241  static int32_t ib555;
242 
243  static int firstCall = 1;
244 
245  l1str *l1rec = l2rec->l1rec;
246 
247  if (firstCall) {
248 
249  /* select 490/555 or 490/565 fit, depending on proximity of */
250  /* sensor bands to fitted bands */
251 
252  if ((ib555 = bindex_get(551)) > 0) {
253  a = 0.1853;
254  b = -1.349;
255  } else if ((ib555 = bindex_get(565)) > 0) {
256  a = 0.1787;
257  b = -1.122;
258  } else {
259  printf("Kd_obpg: incompatible sensor wavelengths (no 555 or 565).\n");
260  exit(1);
261  }
262 
263  if ((ib490 = bindex_get(490)) < 0) {
264  printf("Kd_obpg: incompatible sensor wavelengths (no 490).\n");
265  exit(1);
266  }
267 
268  firstCall = 0;
269  }
270 
271  for (ip = 0; ip < l1rec->npix; ip++) {
272 
273  ipb = l1rec->l1file->nbands*ip;
274 
275  if (l1rec->mask[ip] ||
276  l2rec->nLw[ipb + ib490] <= 0.0 || l2rec->nLw[ipb + ib555] <= 0.0) {
277  k490[ip] = kdbad;
278  l1rec->flags[ip] |= CHLFAIL;
279  l1rec->flags[ip] |= PRODFAIL;
280  } else {
281  k490[ip] = a * pow(l2rec->nLw[ipb + ib490] / l2rec->nLw[ipb + ib555], b);
282  if (k490[ip] > KD_MAX) {
283  k490[ip] = KD_MAX;
284  l1rec->flags[ip] |= PRODWARN;
285  }
286  /* not until reprocessing
287  if (k490[ip] < KD_MIN) {
288  k490[ip] = KD_MIN;
289  l2rec->flags[ip] |= PRODWARN;
290  }
291  */
292  }
293 
294  }
295 }
296 
297 
298 /* ------------------------------------------------------------------- */
299 /* Kd490_morel() - diffuse attenuation at 490 using Morel (2007) */
300 /* */
301 /* Inputs: */
302 /* l2rec - level-2 structure containing one complete scan after */
303 /* atmospheric correction. */
304 /* band - waveband number (0 - nbands-1) at which Kd computed. */
305 /* */
306 /* Outputs: */
307 /* Kd - diffuse attenuation at 490nm, 1 value per pixel. */
308 /* */
309 /* Description: */
310 /* This produces the estimate of diffuse attenation at 490 */
311 /* using the satellite derived chlorophyll. */
312 /* */
313 /* Reference: */
314 /* */
315 /* Morel, A., Y. Huot, B. Gentili, P.J. Werdell, S.B. Hooker (2007). */
316 /* Consistency of products derived from various ocean color sensors: */
317 /* An examination before merging these products and extending their */
318 /* applications, Remote Sensing of Environment, to be submitted. */
319 /* */
320 /* New equation 8 (combined LOV and NOMAD derivation) */
321 /* */
322 /* Original Implementation: B. Franz, February 2007 */
323 
324 /*---------------------------------------------------------------------*/
325 void Kd490_morel(l2str *l2rec, float *Kd) {
326  static float Kw = 0.01660;
327  static float X = 0.077746;
328  static float e = 0.672846;
329 
330  float chl;
331  int32_t ip;
332 
333  l1str *l1rec = l2rec->l1rec;
334 
335  for (ip = 0; ip < l1rec->npix; ip++) {
336 
337  chl = l2rec->chl[ip];
338 
339  if (l1rec->mask[ip] || chl <= 0.0) {
340  Kd[ip] = kdbad;
341  l1rec->flags[ip] |= PRODFAIL;
342  } else {
343  Kd[ip] = Kw + X * pow(chl, e);
344  if (Kd[ip] > KD_MAX) {
345  Kd[ip] = KD_MAX;
346  l1rec->flags[ip] |= PRODWARN;
347  } else
348  if (Kd[ip] < KD_MIN) {
349  Kd[ip] = KD_MIN;
350  l1rec->flags[ip] |= PRODWARN;
351  }
352  }
353  }
354 
355  return;
356 }
357 
358 
359 /* ------------------------------------------------------------------- */
360 /* Kd_PAR_morel() - spectrally integrated attenuation using Morel */
361 /* */
362 /* Inputs: */
363 /* l2rec - level-2 structure containing one complete scan after */
364 /* atmospheric correction. */
365 /* depth - layer depth 1=1/k490, 2=2/k490 */
366 /* */
367 /* Outputs: */
368 /* Kd - Kd(PAR), 1 value per pixel. */
369 /* */
370 /* Description: */
371 /* */
372 /* Reference: */
373 /* */
374 /* Morel, A., Y. Huot, B. Gentili, P.J. Werdell, S.B. Hooker (2007). */
375 /* Consistency of products derived from various ocean color sensors: */
376 /* An examination before merging these products and extending their */
377 /* applications, Remote Sensing of Environment, to be submitted. */
378 /* */
379 /* Original Implementation: B. Franz, November 2006 */
380 
381 /*---------------------------------------------------------------------*/
382 void Kd_PAR_morel(l2str *l2rec, int depth, float *Kd) {
383  static float *Kd490 = NULL;
384 
385  int32_t ip;
386 
387  l1str *l1rec = l2rec->l1rec;
388 
389  if (Kd490 == NULL) {
390  if ((Kd490 = malloc(l1rec->npix * sizeof (float))) == NULL) {
391  printf("-E- %s: Error allocating space for %d records.\n", __FILE__, l1rec->npix);
392  exit(1);
393  }
394  }
395 
396  Kd490_morel(l2rec, Kd490);
397 
398  for (ip = 0; ip < l1rec->npix; ip++) {
399 
400  if (l1rec->mask[ip] || Kd490[ip] <= 0.0) {
401  Kd[ip] = kdbad;
402  l1rec->flags[ip] |= PRODFAIL;
403  } else {
404  switch (depth) {
405  case 1:
406  Kd[ip] = 0.0864 + 0.884 * Kd490[ip] - 0.00137 / Kd490[ip];
407  break;
408  case 2:
409  Kd[ip] = 0.0665 + 0.874 * Kd490[ip] - 0.00121 / Kd490[ip];
410  break;
411  default:
412  printf("-E- %s: Invalid depth for Kd(PAR) (1 or 2).\n", __FILE__);
413  exit(1);
414  break;
415  }
416  }
417  }
418 
419  return;
420 }
421 
422 /* ------------------------------------------------------------------- */
423 /* Zhl_morel() - heated layer depth using Morel (2007) */
424 /* */
425 /* Inputs: */
426 /* l2rec - level-2 structure containing one complete scan after */
427 /* atmospheric correction. */
428 /* */
429 /* Outputs: */
430 /* Zhl - heated layer depth (m) */
431 /* */
432 /* Description: */
433 /* */
434 /* Reference: */
435 /* */
436 /* Morel, A., Y. Huot, B. Gentili, P.J. Werdell, S.B. Hooker (2007). */
437 /* Consistency of products derived from various ocean color sensors: */
438 /* An examination before merging these products and extending their */
439 /* applications, Remote Sensing of Environment, to be submitted. */
440 /* */
441 /* Original Implementation: B. Franz, November 2006 */
442 
443 /*---------------------------------------------------------------------*/
444 void Zhl_morel(l2str *l2rec, float *Zhl) {
445  static float *KdPAR = NULL;
446 
447  int32_t ip;
448  l1str *l1rec = l2rec->l1rec;
449 
450  /* need Kd(PAR) at 2nd optical depth */
451  if (KdPAR == NULL) {
452  if ((KdPAR = malloc(l1rec->npix * sizeof (float))) == NULL) {
453  printf("-E- %s: Error allocating space for %d records.\n", __FILE__, l1rec->npix);
454  exit(1);
455  }
456  }
457  Kd_PAR_morel(l2rec, 2, KdPAR);
458 
459  for (ip = 0; ip < l1rec->npix; ip++) {
460 
461  if (l1rec->mask[ip] || KdPAR[ip] <= 1e-5) {
462  Zhl[ip] = -999;
463  l1rec->flags[ip] |= PRODFAIL;
464  } else
465  Zhl[ip] = 2.0 / KdPAR[ip];
466  }
467 
468  return;
469 }
470 
471 
472 /* ------------------------------------------------------------------- */
473 /* Kd532() - spectral diffuse attenuation using Mueller/Austin&Petzold */
474 /* */
475 /* Inputs: */
476 /* l2rec - level-2 structure containing one complete scan after */
477 /* atmospheric correction. */
478 /* flag - return Kd in meters (1) or 1/meters */
479 /* Kd - diffuse attenuation at 490 nm. */
480 /* Outputs: */
481 /* Kd - diffuse attenuation at 532 nm. */
482 /* */
483 /* Description: */
484 /* This produces the estimate of diffuse attenation at 532 nm from */
485 /* the estimate of diffuse attenuation at 490 nm using the spectral */
486 /* K algorithm of Austin and Petzold. */
487 /* */
488 /* Reference: */
489 /* Austin, R. W and T. J. Petzold, "Spectral Dependence of the */
490 /* Diffuse Attenuation Coefficient of Light in Ocean Waters", */
491 /* SPIE Vol 489 Ocean Optics (1984) pp 168-178 */
492 /* */
493 /* Original Implementation: P. Martinolich, NRL/Stennis, May 2005 */
494 
495 /*---------------------------------------------------------------------*/
496 
497 void Kd532(l2str *l2rec, int flag, float k532[]) {
498  const float M532 = 0.68052;
499  const float KW490 = 0.0224;
500  const float KW532 = 0.05356;
501 
502  static float maxval = 6.4;
503 
504  int ip;
505  float temp;
506 
507  l1str *l1rec = l2rec->l1rec;
508 
509  for (ip = 0; ip < l1rec->npix; ip++) {
510 
511  if (k532[ip] < 0.0)
512  k532[ip] = kdbad;
513  else if (k532[ip] >= maxval)
514  k532[ip] = maxval;
515  else if (l1rec->mask[ip])
516  k532[ip] = kdbad;
517  else {
518  temp = M532 * (k532[ip] - KW490) + KW532;
519  if (flag > 0)
520  k532[ip] = 1.0 / temp;
521  else
522  k532[ip] = temp;
523  if (k532[ip] < 0.0)
524  k532[ip] = kdbad;
525  if (k532[ip] > maxval)
526  k532[ip] = maxval;
527  }
528 
529  }
530 }
531 
559 void Kd_lee(l2str *l2rec, int band, float *Kd)
560 {
561 
562  const float m1 = 4.259;
563  const float m2 = 0.520;
564  const float m3 = -10.800;
565  const float m4 = 0.265;
566 
567  float m0;
568  int ip, ipb;
569 
570  l1str *l1rec = l2rec->l1rec;
571 
572  if (input->iop_opt == IOPNONE) {
573  printf("IOP-based Kd_lee product requires iop model selection (iop_opt). ");
574  printf("Using default model.\n");
575  input->iop_opt = IOPDEFAULT;
576  get_iops(l2rec, input->iop_opt);
577  }
578 
579  for (ip = 0; ip < l1rec->npix; ip++) {
580 
581  ipb = ip * l1rec->l1file->nbands + band;
582 
583  if (l1rec->mask[ip] || l2rec->a[ipb] <= 0.0 || l2rec->bb[ipb] <= 0.0 ) {
584  Kd[ip] = kdbad;
585  l1rec->flags[ip] |= PRODFAIL;
586  } else {
587  m0 = 1.0 + 0.005 * l1rec->solz[ip];
588  Kd[ip] = m0 * l2rec->a[ipb]
589  + m1 * (1.0 - m2 * exp( m3 * l2rec->a[ipb])) * l2rec->bb[ipb]
590  * (1.0 - m4 * (l1rec->sw_bb[ipb] / l2rec->bb[ipb]));
591 
592  if (Kd[ip] > KD_MAX) {
593  Kd[ip] = KD_MAX;
594  l1rec->flags[ip] |= PRODWARN;
595  } else
596  if (Kd[ip] < KD_MIN) {
597  Kd[ip] = KD_MIN;
598  l1rec->flags[ip] |= PRODWARN;
599  }
600  }
601  }
602 
603  return;
604 }
605 
606 /* ------------------------------------------------------------------- */
607 /* Kd_lee_unc() - standard uncertainty for spectral attenuation */
608 /* coefficient per Lee et al (2013) */
609 /* */
610 /* Inputs: */
611 /* l2rec - level-2 structure containing one complete scan after */
612 /* atmospheric correction. */
613 /* band - waveband number (0 - nbands-1) at which nKd computed. */
614 /* */
615 /* Outputs: */
616 /* uKd - standard uncertainty of Kd at specified band, */
617 /* 1 value per pixel. */
618 /* */
619 /* Description: */
620 /* This produces an estimate of the standard uncertainty in the */
621 /* diffuse attenation coffucient computed using Kd_lee algorihtm. */
622 /* It requires the total absorption and backscattering coefficents */
623 /* and their asscociate standard uncertainties as well as the solar */
624 /* zenith angle. */
625 /* */
626 /* Implementation: L. McKinna, November 2023 */
627 /* ------------------------------------------------------------------- */
628 
629 void unc_Kd_lee(l2str *l2rec, int band, float *uKd) {
630 
631  const float m1 = 4.259;
632  const float m2 = 0.520;
633  const float m3 = -10.800;
634  const float m4 = 0.265;
635 
636  float m0;
637  float Kd, dkda, dkdbb;
638  int ip, ipb, iopt_flag;
639 
640  float *a_unc;
641  float *bb_unc;
642 
643  l1str *l1rec = l2rec->l1rec;
644 
645  if (input->iop_opt == IOPNONE) {
646  printf("IOP-based Kd_lee product requires iop model selection (iop_opt). ");
647  printf("Using default model.\n");
648  input->iop_opt = IOPDEFAULT;
649  get_iops(l2rec, input->iop_opt);
650  }
651 
652  switch (input->iop_opt) {
653  case IOPGIOP:
654  iopt_flag = 0;
655  a_unc = giop_get_a_unc_pointer();
656  bb_unc = giop_get_bb_unc_pointer();
657  iopt_flag = 1;
658  break;
659  default:
660  printf("IOP-based uncertainty for Kd_lee requires support iop model (iop_opt). ");
661  printf("Kd uncertanties cannot be estimated, seetting to badval.\n");
662  iopt_flag = 0;
663  break;
664  }
665 
666 
667  for (ip = 0; ip < l1rec->npix; ip++) {
668 
669  ipb = ip * l1rec->l1file->nbands + band;
670 
671  uKd[ip] = BAD_FLT;
672 
673  //Cacluate Kd_lee_unc if iop model produces uncertainties
674  if (iopt_flag) {
675 
676  if (l1rec->mask[ip] || l2rec->a[ipb] <= 0.0 || l2rec->bb[ipb] <= 0.0 ) {
677  uKd[ip] = kdbad;
678  l1rec->flags[ip] |= PRODFAIL;
679 
680  } else {
681  m0 = 1.0 + 0.005 * l1rec->solz[ip];
682  Kd = m0 * l2rec->a[ipb]
683  + m1 * (1.0 - m2 * exp( m3 * l2rec->a[ipb])) * l2rec->bb[ipb]
684  * (1.0 - m4 * (l1rec->sw_bb[ipb] / l2rec->bb[ipb]));
685 
686  dkda = m0 + m1*(1. - m4*(l1rec->sw_bb[ipb] / bb_unc[ipb]))*bb_unc[ipb]
687  * (-m2*m3*exp( m3 * a_unc[ipb]));
688 
689  dkdbb = m1*(1. - m2*exp( m3 * a_unc[ipb]))*(1. - m4*(l1rec->sw_bb[ipb] / l2rec->bb[ipb]));
690 
691  uKd[ip] = sqrt (pow(dkda* a_unc[ipb],2.) + pow(dkdbb*bb_unc[ipb],2));
692 
693  if (Kd > KD_MAX) {
694  uKd[ip] = 0.0;
695  l1rec->flags[ip] |= PRODWARN;
696  } else
697  if (Kd < KD_MIN) {
698  uKd[ip] = 0.0;
699  l1rec->flags[ip] |= PRODWARN;
700  }
701  }
702  }
703 
704  }
705 
706  return;
707 }
708 
709 /* ------------------------------------------------------------------- */
710 /* nKd_lee() - spectral normalized diffuse attenuation using */
711 /* Lin, et. (2016) */
712 /* */
713 /* Inputs: */
714 /* l2rec - level-2 structure containing one complete scan after */
715 /* atmospheric correction. */
716 /* band - waveband number (0 - nbands-1) at which nKd computed. */
717 /* */
718 /* Outputs: */
719 /* nKd - normalized diffuse attenuation at specified band, */
720 /* 1 value per pixel. */
721 /* */
722 /* Description: */
723 /* This produces the estimate of the normalized diffuse attenation at */
724 /* the given band using the satellite derived total absorption and */
725 /* backscattering and solar zenith angle. */
726 /* */
727 /* Reference: */
728 /* Lin, J., Lee, Z.P., Ondrusek, M. and K. Du (2016) "Remote sensing */
729 /* the normalized diffuse attenuation coefficient of downwelling */
730 /* irradiance", J. Geophys. Res. Oceans, 121: 6717-6730 */
731 /* */
732 /* Original Implementation: L. McKinna, June 2017 */
733 
734 /*---------------------------------------------------------------------*/
735 void nKd_lin(l2str *l2rec, int band, float *nKd) {
736 
737  const float m1 = 4.18;
738  const float m2 = 0.52;
739  const float m3 = -10.80;
740 
741  const float n1 = 7.41;
742  const float n2 = 0.13;
743  const float n3 = 4.36;
744  const float n4 = -0.0051;
745  const float n5 = 0.0094;
746  const float n6 = 1.78;
747  const float p = 1.14;
748  const float radCon = M_PI / 180.0; //degrees -> radians conversion
749 
750  float m0, sasw, D0, f, lambda, Kd;
751 
752  int ip, ipb;
753 
754  l1str *l1rec = l2rec->l1rec;
755 
756  if (input->iop_opt == IOPNONE) {
757  printf("IOP-based Kd_lee product requires iop model selection (iop_opt). ");
758  printf("Using default model.\n");
759  input->iop_opt = IOPDEFAULT;
760  get_iops(l2rec, input->iop_opt);
761  }
762 
763  for (ip = 0; ip < l1rec->npix; ip++) {
764 
765  ipb = ip * l1rec->l1file->nbands + band;
766  lambda = l2rec->l1rec->l1file->fwave[band];
767  sasw = asin(sin(radCon * l1rec->solz[ip]) / 1.34);
768  f = n1 * (n2 - exp(-n3 * (lambda / 490.))) - (n4 * (lambda / 490.) + n5) * exp(n6 * (l1rec->solz[ip] / 30.));
769  D0 = (f / cos(sasw)) + p * (1.0 - f);
770 
771  if (l1rec->mask[ip] ||
772  l2rec->a[ipb] <= 0.0 || l2rec->bb[ipb] <= 0.0) {
773  nKd[ip] = kdbad;
774  l1rec->flags[ip] |= PRODFAIL;
775  } else {
776  m0 = 1.0 + 0.005 * l1rec->solz[ip];
777  Kd = m0 * l2rec->a[ipb]
778  + m1 * (1.0 - m2 * exp(m3 * l2rec->a[ipb])) * l2rec->bb[ipb];
779 
780  nKd[ip] = Kd / D0;
781 
782  if (nKd[ip] > KD_MAX) {
783  nKd[ip] = KD_MAX;
784  l1rec->flags[ip] |= PRODWARN;
785  } else
786  if (nKd[ip] < KD_MIN) {
787  nKd[ip] = KD_MIN;
788  l1rec->flags[ip] |= PRODWARN;
789  }
790  }
791  }
792 
793  return;
794 }
795 
796 /* ------------------------------------------------------------------- */
797 /* Kd_PAR_lee() - spectrally integrated attenuation using ZP Lee */
798 /* */
799 /* Inputs: */
800 /* l2rec - level-2 structure containing one complete scan after */
801 /* atmospheric correction. */
802 /* */
803 /* Outputs: */
804 /* Kd - Kd(PAR), 1 value per pixel. */
805 /* */
806 /* Description: */
807 /* */
808 /* Reference: */
809 /* */
810 /* Original Implementation: B. Franz, September 2007 */
811 
812 /*---------------------------------------------------------------------*/
813 void Kd_PAR_lee(l2str *l2rec, float *Kd) {
814  void Zphotic_lee(l2str *l2rec, l2prodstr *p, float Zp[]);
815 
816  float *Zp = Kd;
817  l2prodstr p;
818  int32_t ip;
819 
820  // get first optical depth from Lee
821  p.cat_ix = CAT_Zphotic_lee;
822  p.prod_ix = -3;
823  Zphotic_lee(l2rec, &p, Zp);
824 
825  // invert to get Kd(PAR) at 1st optical depth
826  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
827  if (Zp[ip] > 0.0)
828  Kd[ip] = 1.0 / Zp[ip];
829  else {
830  Kd[ip] = kdbad;
831  l2rec->l1rec->flags[ip] |= PRODFAIL;
832  }
833  }
834 }
835 
836 /* ------------------------------------------------------------------- */
837 /* Kd_jamet() - spectral diffuse attenuation using Jamet et al, (2012) */
838 /* */
839 /* Inputs: */
840 /* l2rec - level-2 structure containing one complete scan after */
841 /* atmospheric correction. */
842 /* band - waveband number (0 - nbands-1) at which Kd computed. */
843 /* */
844 /* Outputs: */
845 /* Kd - diffuse attenuation at specified band, 1 value per pixel. */
846 /* */
847 /* Description: */
848 /* This derives the Kd using the full set of Rrs available and only */
849 /* for the SeaWiFS instrument (at this time). A neural network is */
850 /* used to perform this. */
851 /* */
852 /* Reference: */
853 /* Jamet, C., H. Loisel, and D. Dessailly (2012), Retrieval of the */
854 /* spectral diffuse attenuation coefficient Kd(lambda) in open and */
855 /* coastal ocean waters using a neural network inversion, J Geophys. */
856 /* Res., 117, C10023, doi: 10.1029/2012JC008076. */
857 /* */
858 /* Original Implementation: C. Jamet */
859 
860 /*---------------------------------------------------------------------*/
861 void Kd_jamet(l2str *l2rec, int band, float *Kd) {
862  static int firstCall = 1;
863  int ip, ipb;
864  static float alg_lambda[16], alg_inp[16];
865  int32 i, j, bad_prd;
866  static float *b1, *w1, *w2, *mean, *sd, b2;
867  static float *xnorm, *alayer;
868  static int32 nrrs, n_input, n_norm, n_nodes;
869  static int32 *bnd_ptr;
870  float y;
871 
872  filehandle *l1file = l2rec->l1rec->l1file;
873 
874  /*
875  * Initialize the coefficients for the computation
876  */
877  if (firstCall) {
878  FILE *fp;
879  char *filedir;
880  char filename[FILENAME_MAX];
881  char line [801];
882  int com_ln;
883  int poub;
884 
885  firstCall = 0;
886  if ((bnd_ptr = (int32 *) calloc(l1file->nbands, sizeof (int32))) == NULL) {
887  printf("-E- : Error allocating memory to bnd_ptr in get_Kd\n");
888  exit(FATAL_ERROR);
889  }
890  if ((filedir = getenv("OCDATAROOT")) == NULL) {
891  printf("-E- %s, %d: OCDATAROOT env variable undefined.\n",
892  __FILE__, __LINE__);
893  exit(1);
894  }
895  strcpy(filename, filedir);
896  switch (l1file->sensorID) {
897  case SEAWIFS:
898  strcat(filename, "/seawifs/iop/kd_jamet/seawifs_kd_jamet.dat");
899  break;
900  case MERIS:
901  strcat(filename, "/meris/iop/kd_jamet/meris_kd_jamet.dat");
902  break;
903  default:
904  printf(
905  "-E- %s, %d: Kd_jamet can only be generated for the SEAWIFS instrument\n",
906  __FILE__, __LINE__);
907  exit(1);
908  break;
909  }
910  printf("Loading Kd coefficient table %s\n", filename);
911 
912  if ((fp = fopen(filename, "r")) == NULL) {
913  fprintf(stderr, "-E- %s, %d: unable to open %s for reading\n",
914  __FILE__, __LINE__, filename);
915  exit(1);
916  }
917  /*
918  * read the # bands for the Rrs and their values
919  */
920  com_ln = 1;
921  while (com_ln == 1) {
922  if (fgets(line, 800, fp) == NULL) {
923  fprintf(stderr, "-E- %s, %d: Error reading %s\n",
924  __FILE__, __LINE__, filename);
925  exit(1);
926  }
927  if (line[0] != ';') com_ln = 0;
928  }
929  if (sscanf(line, "%d", &nrrs) != 1) {
930  fprintf(stderr, "-E- %s, %d: Error reading %s\n",
931  __FILE__, __LINE__, filename);
932  exit(1);
933  }
934  for (i = 0; i < nrrs; i++) {
935  fgets(line, 800, fp);
936  sscanf(line, "%f", &alg_lambda[i]);
937  }
938  n_input = nrrs + 1;
939  n_norm = n_input + 1;
940  /*
941  * get the # nodes in the hidden layer
942  */
943  com_ln = 1;
944  while (com_ln == 1) {
945  fgets(line, 800, fp);
946  if (line[0] != ';') com_ln = 0;
947  }
948  sscanf(line, "%d", &n_nodes);
949 
950  b1 = malloc(n_nodes * sizeof ( float));
951  w1 = malloc(n_nodes * n_norm * sizeof ( float));
952  w2 = malloc(n_nodes * sizeof ( float));
953  mean = malloc(n_norm * sizeof ( float));
954  sd = malloc(n_norm * sizeof ( float));
955  xnorm = malloc(n_input * sizeof ( float));
956  alayer = malloc(n_nodes * sizeof ( float));
957 
958  if ((b1 == NULL) || (w1 == NULL) || (w2 == NULL) ||
959  (mean == NULL) || (sd == NULL) ||
960  (xnorm == NULL) || (alayer == NULL)) {
961  fprintf(stderr, "-E- %s, %d: allocation of Kd weights failed\n",
962  __FILE__, __LINE__);
963  exit(1);
964  }
965  /*
966  * get the weights for the hidden and output layers
967  */
968  com_ln = 1;
969  while (com_ln == 1) {
970  fgets(line, 800, fp);
971  if (line[0] != ';') com_ln = 0;
972  }
973  /* ( this line is a throw-away line) */
974  for (i = 0; i < n_nodes; i++) {
975  fgets(line, 800, fp);
976  sscanf(line, "%d %d %f", &poub, &poub, &b1[i]);
977  }
978  fgets(line, 800, fp);
979  sscanf(line, "%d %d %f", &poub, &poub, &b2);
980  for (j = 0; j < n_input; j++)
981  for (i = 0; i < n_nodes; i++) {
982  fgets(line, 800, fp);
983  sscanf(line, "%d %d %f", &poub, &poub, (w1 + i + n_nodes * j));
984  }
985  for (j = 0; j < n_nodes; j++) {
986  fgets(line, 800, fp);
987  sscanf(line, "%d %d %f", &poub, &poub, &w2[j]);
988  }
989  /*
990  * lastly, read the mean and standard deviation data
991  */
992  com_ln = 1;
993  while (com_ln == 1) {
994  fgets(line, 800, fp);
995  if (line[0] != ';') com_ln = 0;
996  }
997  sscanf(line, "%f", &mean[0]);
998  for (i = 1; i < n_norm; i++) {
999  fgets(line, 800, fp);
1000  sscanf(line, "%f", &mean[i]);
1001  }
1002  com_ln = 1;
1003  while (com_ln == 1) {
1004  fgets(line, 800, fp);
1005  if (line[0] != ';') com_ln = 0;
1006  }
1007  sscanf(line, "%f", &sd[0]);
1008  for (i = 1; i < n_norm; i++) {
1009  fgets(line, 800, fp);
1010  sscanf(line, "%f", &sd[i]);
1011  }
1012  fclose(fp);
1013  /*
1014  * set up correct band locations
1015  */
1016  printf("Kd_jamet band assignment:\n");
1017  printf("Alg wave inst wave inst band\n");
1018  for (i = 0; i < nrrs; i++) {
1019  bnd_ptr[i] = windex(alg_lambda[i], l1file->fwave, l1file->nbands);
1020  printf("%7.2f %7.2f %3d\n", alg_lambda[i],
1021  l1file->fwave[bnd_ptr[i]], bnd_ptr[i]);
1022  }
1023  printf("\n");
1024  }
1025  /*
1026  * derive the Kd
1027  */
1028  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
1029  ipb = ip * l1file->nbands;
1030  /*
1031  * Normalize the Rrs, lambda
1032  */
1033  bad_prd = 0;
1034  if (l2rec->l1rec->mask[ip]) bad_prd = 1;
1035  /*
1036  for( i = 0; i < nrrs; i++ )
1037  {
1038  alg_inp[i] = l2rec->Rrs[ ipb + bnd_ptr[i] ];
1039  if( i <= 1 ) {
1040  if( alg_inp[i] <= -0.01 ) bad_prd = 1;
1041  }
1042  else {
1043  if( alg_inp[i] <= 0.0 ) bad_prd = 1;
1044  }
1045  }
1046  */
1047  for (i = 0; i < nrrs; i++) {
1048  alg_inp[i] = l2rec->Rrs[ ipb + bnd_ptr[i] ];
1049  if (alg_inp[i] < BAD_FLT + 1) bad_prd = 1;
1050  }
1051  alg_inp[nrrs] = l1file->fwave[band];
1052 
1053  if (bad_prd == 1) {
1054  Kd[ip] = kdbad;
1055  l2rec->l1rec->flags[ip] |= PRODFAIL;
1056  } else {
1057  for (i = 0; i < n_input; i++)
1058  xnorm[i] = ((2. / 3.) *(alg_inp[i] - mean[i])) / sd[i];
1059  /*
1060  * Apply the nn algorithm
1061  */
1062  for (i = 0; i < n_nodes; i++) {
1063  alayer[i] = 0.0;
1064  for (j = 0; j < n_input; j++) {
1065  alayer[i] += (xnorm[j] * *(w1 + i + n_nodes * j));
1066  }
1067  alayer[i] = 1.715905 * (float) tanh((2. / 3.) *
1068  (double) (alayer[i] + b1[i]));
1069  }
1070 
1071  for (j = 0, y = 0.; j < n_nodes; j++)
1072  y += (alayer[j] * w2[j]);
1073  /*
1074  * De-normalize the log(Kd) and make Kd
1075  */
1076  y = 1.5 * (y + b2) * sd[n_input] + mean[n_input];
1077  *(Kd + ip) = (float) pow(10., (double) y);
1078  }
1079  }
1080  return;
1081 }
1082 
1083 /* ------------------------------------------------------------------- */
1084 /* Kd_rhos() - spectral diffuse attenuation using Stumpf, et. (2013) */
1085 /* */
1086 /* Inputs: */
1087 /* l2rec - level-2 structure containing one complete scan after */
1088 /* atmospheric correction. */
1089 /* band - waveband number (0 - nbands-1) at which Kd computed. */
1090 /* */
1091 /* Outputs: */
1092 /* Kd - diffuse attenuation at specified band, 1 value per pixel. */
1093 /* */
1094 /* Description: */
1095 /* This produces the estimate of diffuse attenation at the given band */
1096 /* using the satellite derived total absorption and backscattering. */
1097 /* */
1098 /* Reference: */
1099 /* Original Implementation: R. Healy, NASA-GSFC, August 2015 */
1100 
1101 /*---------------------------------------------------------------------*/
1102 void Kd_rhos(l2str *l2rec, float *Kd) {
1103 
1104  const float factor = 0.7;
1105 
1106  static int ib443, ib490, ib620, ib665, ib645, ib469, ib859, ib865, firstCall = 1;
1107  int ip, ipb;
1108 
1109  if (firstCall == 1) {
1110 
1111  ib645 = bindex_get(645);
1112  ib469 = bindex_get(469);
1113  ib859 = bindex_get(859);
1114  firstCall = 0;
1115 
1116  ib443 = -1;
1117  if (ib645 < 0 || ib469 < 0 || ib859 < 0) {
1118  // try the Meris wavelengths
1119 
1120  ib443 = bindex_get(443);
1121  ib490 = bindex_get(490);
1122  ib620 = bindex_get(620);
1123  ib665 = bindex_get(665);
1124  ib865 = bindex_get(865);
1125  firstCall = 0;
1126 
1127 
1128  if (ib443 < 0 || ib490 < 0 || ib620 < 0 || ib665 < 0 || ib865 < 0) {
1129  printf("Kd_rhos: incompatible sensor wavelengths for this algorithm\n");
1130  exit(1);
1131  } else {
1132  printf("Kd_rhos: Using the Meris band averaging for this product\n");
1133  }
1134 
1135  }
1136 
1137  }
1138 
1139  if (ib443 < 0) {
1140  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
1141 
1142  ipb = l2rec->l1rec->l1file->nbands*ip;
1143 
1144  // if (l2rec->mask[ip] ||
1145  if (l2rec->Rrs[ipb + ib645] <= 0.0 || l2rec->Rrs[ipb + ib469] <= 0.0 || l2rec->Rrs[ipb + ib859] <= 0.0) {
1146  Kd[ip] = kdbad;
1147  l2rec->l1rec->flags[ip] |= PRODFAIL;
1148  } else {
1149  Kd[ip] = factor * (l2rec->l1rec->rhos[ipb + ib645] - l2rec->l1rec->rhos[ipb + ib859])
1150  / (l2rec->l1rec->rhos[ipb + ib469] - l2rec->l1rec->rhos[ipb + ib859]);
1151  }
1152  }
1153  } else {
1154  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
1155 
1156  ipb = l2rec->l1rec->l1file->nbands*ip;
1157 
1158  // if (l2rec->mask[ip] ||
1159  if (l2rec->Rrs[ipb + ib665] <= 0.0 || l2rec->Rrs[ipb + ib620] <= 0.0 || l2rec->Rrs[ipb + ib443] <= 0.0 || l2rec->Rrs[ipb + ib490] <= 0.0 || l2rec->Rrs[ipb + ib865] <= 0.0) {
1160  Kd[ip] = kdbad;
1161  l2rec->l1rec->flags[ip] |= PRODFAIL;
1162  } else {
1163  Kd[ip] = factor * ((l2rec->l1rec->rhos[ipb + ib620] + l2rec->l1rec->rhos[ipb + ib665]) / 2 - l2rec->l1rec->rhos[ipb + ib865])
1164  / ((l2rec->l1rec->rhos[ipb + ib443] + l2rec->l1rec->rhos[ipb + ib490]) / 2 - l2rec->l1rec->rhos[ipb + ib865]);
1165  }
1166 
1167  }
1168  }
1169  return;
1170 }
1171 
1172 /* ------------------------------------------------------------------- */
1173 /* get_spectral_Kd() - untilty function to generate spectral Kd pointer*/
1174 /* */
1175 /* Inputs: */
1176 /* l2rec - level-2 structure containing one complete scan after */
1177 /* atmospheric correction. */
1178 /* l2prodstr - level-2 product structure containing information */
1179 /* about requested product(s) */
1180 /* */
1181 /* Outputs: */
1182 /* Kd - diffuse attenuation at specified band, 1 value per pixel. */
1183 /* */
1184 /*---------------------------------------------------------------------*/
1185 
1186 void get_spectral_Kd(l2str *l2rec,l2prodstr *p) {
1187 
1188  l1str *l1rec = l2rec->l1rec;
1189 
1190  int32 iw,ip;
1191  int32 npix = l1rec->npix;
1192  int32 nwave = l1rec->l1file->nbands;
1193 
1194  float temp;
1195  float *Kd_temp;
1196  Kd_temp = (float*) calloc(npix, sizeof (float));
1197 
1198  Kd_spectral = (float*)calloc(npix*nwave,sizeof (float));
1199 
1200  for (iw=0;iw<nwave; iw++){
1201 
1202  //Switch between requested product
1203  switch (p->cat_ix) {
1204  case CAT_Kd_lee:
1205  Kd_lee(l2rec, iw, Kd_temp);
1206  break;
1207  case CAT_Kd_unc_lee:
1208  unc_Kd_lee(l2rec, iw, Kd_temp);
1209  break;
1210  default:
1211  printf("Error: %s : 3d wavelength Kd not supported for model: %d\n", __FILE__, p->cat_ix);
1212  exit(1);
1213  break;
1214  }
1215 
1216  for (ip=0;ip<npix; ip++){
1217  temp = Kd_temp[ip];
1218  Kd_spectral[ip*nwave +iw] = temp;
1219  }
1220 
1221  }
1222 
1223  free(Kd_temp);
1224 }
1225 
1226 /*----------------------------------------------------------------------*/
1227 /* extract_band_3d() - utility function to extract selected band subset */
1228 
1229 /*----------------------------------------------------------------------*/
1230 static void extract_band_3d(float *in_buf, float *out_buf, int numPixels, int numBands) {
1231  float * out_ptr = out_buf;
1232  for (int pix = 0; pix < numPixels; pix++) {
1233  float *in_ptr = in_buf + pix * numBands;
1234  for (int band_3d = 0; band_3d < input->nwavelengths_3d; band_3d++) {
1235  int band = input->wavelength_3d_index[band_3d];
1236  *out_ptr = in_ptr[band];
1237  out_ptr++;
1238  }
1239  }
1240 }
1241 
1242 /* ------------------------------------------------------------------- */
1243 /* get_Kd() - l2_hdf_generic interface for Kd */
1244 
1245 /* ------------------------------------------------------------------- */
1246 void get_Kd(l2str *l2rec, l2prodstr *p, float prod[]) {
1247 
1248 
1249  if(p->rank == 3) {
1250 
1251  get_spectral_Kd(l2rec,p);
1252 
1253  switch (p->cat_ix) {
1254 
1255  case CAT_Kd_lee:
1256  extract_band_3d(Kd_spectral, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
1257  break;
1258  case CAT_Kd_unc_lee:
1259  extract_band_3d(Kd_spectral, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
1260  break;
1261  default:
1262  printf("Error: %s : 3d wavelenght extract not support for Kd model: %d\n", __FILE__, p->cat_ix);
1263  exit(1);
1264  break;
1265  }
1266 
1267  } else {
1268 
1269  switch (p->cat_ix) {
1270  case CAT_Kd_mueller:
1271  Kd490_mueller(l2rec, prod);
1272  break;
1273  case CAT_Kd_obpg:
1274  Kd490_obpg(l2rec, prod);
1275  break;
1276  case CAT_Kd_KD2:
1277  Kd490_KD2(l2rec, prod);
1278  break;
1279  case CAT_Kd_unc_KD2:
1280  unc_Kd490_KD2(l2rec, prod);
1281  break;
1282  case CAT_Kd_lee:
1283  Kd_lee(l2rec, p->prod_ix, prod);
1284  break;
1285  case CAT_Kd_unc_lee:
1286  unc_Kd_lee(l2rec, p->prod_ix, prod);
1287  break;
1288  case CAT_Kd_532:
1289  Kd490_mueller(l2rec, prod);
1290  Kd532(l2rec, p->prod_ix, prod);
1291  break;
1292  case CAT_Kd_morel:
1293  Kd490_morel(l2rec, prod);
1294  break;
1295  case CAT_Kd_jamet:
1296  Kd_jamet(l2rec, p->prod_ix, prod);
1297  break;
1298  case CAT_Kd_rhos:
1299  Kd_rhos(l2rec, prod);
1300  break;
1301  case CAT_KPAR_morel:
1302  Kd_PAR_morel(l2rec, p->prod_ix, prod);
1303  break;
1304  case CAT_KPAR_lee:
1305  Kd_PAR_lee(l2rec, prod);
1306  break;
1307  case CAT_Zhl_morel:
1308  Zhl_morel(l2rec, prod);
1309  break;
1310  case CAT_nKd_lin:
1311  nKd_lin(l2rec, p->prod_ix, prod);
1312  break;
1313  default:
1314  printf("Error: %s : Unknown product specifier: %d\n", __FILE__, p->cat_ix);
1315  exit(1);
1316  break;
1317  }
1318  }
1319 
1320  return;
1321 }
1322 
1323 void Kd490_unc(l2str *l2rec,int32_t ip) {
1324  static int32_t *w = NULL;
1325  static float *a = NULL;
1326  static int ib1 = -1;
1327  static int ib2 = -1;
1328  float fit_unc=0.1;
1329 
1330  int32_t ipb;
1331  float R;
1332 
1333  l1str *l1rec = l2rec->l1rec;
1334 
1335  uncertainty_t *uncertainty=l2rec->l1rec->uncertainty;
1336  float *dRrs, *dkd490, *covariance_matrix, kd_temp,tmpderiv,*Rrs;
1337  static int nbands=0;
1338 
1339  if(uncertainty){
1340  dRrs=uncertainty->dRrs;
1341  dkd490=&uncertainty->dkd490;
1342  if(input->proc_uncertainty==2)
1343  covariance_matrix=uncertainty->covaraince_matrix;
1344  else
1345  covariance_matrix=uncertainty->pixel_covariance;
1346  }
1347 
1348  if (w == NULL) {
1349  nbands=l2rec->l1rec->l1file->nbands;
1350  w = input->kd2w;
1351  a = input->kd2c;
1352  if (w[0] < 0 || w[1] < 0) {
1353  printf("Kd490_KD2: algorithm coefficients not provided for this sensor.\n");
1354  exit(1);
1355  }
1356  ib1 = bindex_get(w[0]);
1357  ib2 = bindex_get(w[1]);
1358  if (ib1 < 0 || ib2 < 0) {
1359  printf("Kd490_KD2: incompatible sensor wavelengths for this algorithm\n");
1360  exit(1);
1361  }
1362  }
1363 
1364  ipb = l1rec->l1file->nbands*ip;
1365  Rrs=&l2rec->Rrs[ipb];
1366 
1367  if (l1rec->mask[ip] ||
1368  l2rec->Rrs[ipb + ib1] <= 0.0 || l2rec->Rrs[ipb + ib2] <= 0.0) {
1369  *dkd490 = 0;
1370  } else {
1371  R = log10(l2rec->Rrs[ipb + ib1] / l2rec->Rrs[ipb + ib2]);
1372  if (isnan(R)) {
1373  *dkd490 = 0;
1374  } else {
1375  kd_temp = a[0] + pow(10.0, a[1] + R * (a[2] + R * (a[3] + R * (a[4] + R * a[5]))));
1376 
1377  if (kd_temp > KD_MAX) {
1378  *dkd490 = 0;
1379  }
1380  else{
1381  tmpderiv= (kd_temp-a[0])* (a[2]+2*a[3]*R+3*a[4]*R*R+4*a[5]*R*R*R)*log(10);
1382  tmpderiv/=(log(10)*Rrs[ib1]/Rrs[ib2]);
1383 
1384  *dkd490=pow(dRrs[ib1]/Rrs[ib2],2)+ pow(dRrs[ib2]*Rrs[ib1]/(Rrs[ib2]*Rrs[ib2]),2);
1385 
1386  *dkd490-=2*Rrs[ib1]/(Rrs[ib2]*Rrs[ib2]*Rrs[ib2])*covariance_matrix[ib1*nbands+ib2];
1387 
1388  *dkd490=fabs(tmpderiv)* sqrt(*dkd490);
1389 
1390  *dkd490=sqrt(*dkd490**dkd490+fit_unc*fit_unc*kd_temp*kd_temp)/kd_temp;
1391  }
1392  }
1393  }
1394 
1395 }
1396 
1397 /* =================================================================== */
1398 /* defunct versions */
1399 /* =================================================================== */
1400 
1401 
1402 /* ------------------------------------------------------------------- */
1403 /* Kd_morel() - spectral diffuse attenuation using Morel (2007) */
1404 /* */
1405 /* Inputs: */
1406 /* l2rec - level-2 structure containing one complete scan after */
1407 /* atmospheric correction. */
1408 /* band - waveband number (0 - nbands-1) at which Kd computed. */
1409 /* */
1410 /* Outputs: */
1411 /* Kd - diffuse attenuation at specified band, 1 value per pixel. */
1412 /* */
1413 /* Description: */
1414 /* This produces the estimate of diffuse attenation at the given band */
1415 /* using the satellite derived chlorophyll. */
1416 /* */
1417 /* Reference: */
1418 /* */
1419 /* Morel, A., Y. Huot, B. Gentili, P.J. Werdell, S.B. Hooker (2007). */
1420 /* Consistency of products derived from various ocean color sensors: */
1421 /* An examination before merging these products and extending their */
1422 /* applications, Remote Sensing of Environment, to be submitted. */
1423 /* */
1424 /* Equation 8 */
1425 /* */
1426 /* Original Implementation: B. Franz, August 2006 */
1427 /*---------------------------------------------------------------------*/
1428 /*
1429 void Kd_morel(l2str *l2rec, int band, float *Kd)
1430 {
1431  typedef struct kd_morel_table {
1432  float wave;
1433  float Kw;
1434  float X;
1435  float e;
1436  } kdtabstr;
1437 
1438  static kdtabstr *kdtab = NULL;
1439  static int ntab = 0;
1440 
1441  float chl;
1442  float Kd1, Kd2;
1443  float wt;
1444  float wave;
1445  int32_t ip, itab, i1, i2;
1446 
1447  if (kdtab == NULL) {
1448  FILE *fp = NULL;
1449  char *tmp_str;
1450  char file [FILENAME_MAX] = "";
1451  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
1452  printf("OCDATAROOT environment variable is not defined.\n");
1453  exit(1);
1454  }
1455  strcpy(file,tmp_str); strcat(file,"/common/kd_morel.dat");
1456  if ( (fp = fopen(file,"r")) == NULL ) {
1457  printf("-E- %s: Error opening file %s.\n",__FILE__,file);
1458  exit(1);
1459  }
1460  printf("\nLoading Morel spectral Kd table\n %s\n",file);
1461  fscanf(fp,"%d",&ntab);
1462  if (ntab <= 0) {
1463  printf("-E- %s: Error reading %s.\n",__FILE__,file);
1464  exit(1);
1465  }
1466  if ((kdtab = malloc(ntab*sizeof(kdtabstr))) == NULL) {
1467  printf("-E- %s: Error allocating space for %d records.\n",__FILE__,ntab);
1468  exit(1);
1469  }
1470  for (itab=0; itab<ntab; itab++) {
1471  fscanf(fp,"%f %f %f %f",
1472  &kdtab[itab].wave,
1473  &kdtab[itab].Kw,
1474  &kdtab[itab].X,
1475  &kdtab[itab].e);
1476  printf("%f %f %f %f\n",
1477  kdtab[itab].wave,
1478  kdtab[itab].Kw,
1479  kdtab[itab].X,
1480  kdtab[itab].e);
1481  }
1482  fclose(fp);
1483  printf("\n");
1484  }
1485 
1486  if (band >= 0)
1487  wave = l2rec->fwave[band];
1488  else
1489  wave = 490.;
1490 
1491  for (itab=0; itab<ntab; itab++)
1492  if (wave <= kdtab[itab].wave)
1493  break;
1494 
1495  if (wave == kdtab[itab].wave) {
1496  i1 = itab;
1497  i2 = itab;
1498  } else if (itab == ntab) {
1499  i1 = itab-2;
1500  i2 = itab-1;
1501  wt = (wave-kdtab[i1].wave)/(kdtab[i2].wave - kdtab[i1].wave);
1502  } else {
1503  i1 = itab;
1504  i2 = itab+1;
1505  wt = (wave-kdtab[i1].wave)/(kdtab[i2].wave - kdtab[i1].wave);
1506  }
1507 
1508  for (ip=0; ip<l2rec->npix; ip++) {
1509 
1510  chl = l2rec->chl[ip];
1511 
1512  if (l2rec->mask[ip] || chl <= 0.0) {
1513  Kd[ip] = kdbad;
1514  l2rec->flags[ip] |= PRODFAIL;
1515  } else {
1516  if (i1 == i2) {
1517  Kd[ip] = kdtab[i1].Kw + kdtab[i1].X * pow(chl,kdtab[i1].e);
1518  } else {
1519  Kd1 = kdtab[i1].Kw + kdtab[i1].X * pow(chl,kdtab[i1].e);
1520  Kd2 = kdtab[i2].Kw + kdtab[i2].X * pow(chl,kdtab[i2].e);
1521  Kd[ip] = Kd1 + (Kd2-Kd1)*wt;
1522  }
1523  if (Kd[ip] > KD_MAX) {
1524  Kd[ip] = KD_MAX;
1525  l2rec->flags[ip] |= PRODWARN;
1526  } else
1527  if (Kd[ip] < KD_MIN) {
1528  Kd[ip] = KD_MIN;
1529  l2rec->flags[ip] |= PRODWARN;
1530  }
1531  }
1532  }
1533 
1534  return;
1535 }
1536  */
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
j
int j
Definition: decode_rs.h:73
get_spectral_Kd
void get_spectral_Kd(l2str *l2rec, l2prodstr *p)
Definition: get_Kd.c:1186
get_Kd
void get_Kd(l2str *l2rec, l2prodstr *p, float prod[])
Definition: get_Kd.c:1246
CAT_KPAR_lee
#define CAT_KPAR_lee
Definition: l2prod.h:164
mean
float mean(float *xs, int sample_size)
Definition: numerical.c:81
extract_band_3d
float * extract_band_3d(float *in_buf, float *out_buf)
Definition: prodgen.c:67
CHLFAIL
#define CHLFAIL
Definition: l2_flags.h:26
NULL
#define NULL
Definition: decode_rs.h:63
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
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
Kd532
void Kd532(l2str *l2rec, int flag, float k532[])
Definition: get_Kd.c:497
PRODWARN
#define PRODWARN
Definition: l2_flags.h:13
MERIS
#define MERIS
Definition: sensorDefs.h:22
Kd_PAR_lee
void Kd_PAR_lee(l2str *l2rec, float *Kd)
Definition: get_Kd.c:813
IOPDEFAULT
#define IOPDEFAULT
Definition: l12_parms.h:76
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:41
lambda
data_t lambda[NROOTS+1]
Definition: decode_rs.h:75
Kd_PAR_morel
void Kd_PAR_morel(l2str *l2rec, int depth, float *Kd)
Definition: get_Kd.c:382
IOPNONE
#define IOPNONE
Definition: l12_parms.h:67
KD_MIN
#define KD_MIN
Definition: get_Kd.c:12
input
instr * input
Definition: main_l2binmatch.cpp:27
M_PI
#define M_PI
Definition: dtranbrdf.cpp:19
names::if
character(len=1000) if
Definition: names.f90:13
IOPGIOP
#define IOPGIOP
Definition: l12_parms.h:74
Kd490_morel
void Kd490_morel(l2str *l2rec, float *Kd)
Definition: get_Kd.c:325
l12_proto.h
KD_MAX
#define KD_MAX
Definition: get_Kd.c:11
giop_get_a_unc_pointer
float * giop_get_a_unc_pointer()
Definition: giop.c:3814
bindex_get
int bindex_get(int32_t wave)
Definition: windex.c:45
Kd490_mueller
void Kd490_mueller(l2str *l2rec, float k490[])
Definition: get_Kd.c:184
f
double precision function f(R1)
Definition: tmd.lp.f:1454
Kd490_obpg
void Kd490_obpg(l2str *l2rec, float k490[])
Definition: get_Kd.c:235
CAT_Kd_unc_lee
#define CAT_Kd_unc_lee
Definition: l2prod.h:547
Kd_jamet
void Kd_jamet(l2str *l2rec, int band, float *Kd)
Definition: get_Kd.c:861
CAT_KPAR_morel
#define CAT_KPAR_morel
Definition: l2prod.h:144
color_dtdb.y
y
Definition: color_dtdb.py:430
nKd_lin
void nKd_lin(l2str *l2rec, int band, float *nKd)
Definition: get_Kd.c:735
get_iops
void get_iops(l2str *l2rec, int32_t iop_opt)
Definition: convl12.c:118
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
CAT_Kd_morel
#define CAT_Kd_morel
Definition: l2prod.h:142
create_images.Rrs
def Rrs
Definition: create_images.py:92
unc_Kd490_KD2
void unc_Kd490_KD2(l2str *l2rec, float *uKd)
Definition: get_Kd.c:95
filename
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
fix_mac_rpath.line
line
Definition: fix_mac_rpath.py:56
CAT_Kd_mueller
#define CAT_Kd_mueller
Definition: l2prod.h:38
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
b
data_t b[NROOTS+1]
Definition: decode_rs.h:77
CAT_Kd_KD2
#define CAT_Kd_KD2
Definition: l2prod.h:150
CAT_Kd_jamet
#define CAT_Kd_jamet
Definition: l2prod.h:283
CAT_Kd_532
#define CAT_Kd_532
Definition: l2prod.h:163
CAT_Zphotic_lee
#define CAT_Zphotic_lee
Definition: l2prod.h:160
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
OCTS
#define OCTS
Definition: sensorDefs.h:14
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
fabs
#define fabs(a)
Definition: misc.h:93
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
giop_get_bb_unc_pointer
float * giop_get_bb_unc_pointer()
Definition: giop.c:3819
Zphotic_lee
void Zphotic_lee(l2str *l2rec, l2prodstr *p, float Zp[])
Definition: photic_depth.c:37
CAT_Kd_obpg
#define CAT_Kd_obpg
Definition: l2prod.h:109
Kd_rhos
void Kd_rhos(l2str *l2rec, float *Kd)
Definition: get_Kd.c:1102
R
#define R
Definition: make_L3_v1.1.c:96
CAT_Kd_unc_KD2
#define CAT_Kd_unc_KD2
Definition: l2prod.h:543
CAT_Zhl_morel
#define CAT_Zhl_morel
Definition: l2prod.h:145
unc_Kd_lee
void unc_Kd_lee(l2str *l2rec, int band, float *uKd)
Definition: get_Kd.c:629
SEAWIFS
#define SEAWIFS
Definition: sensorDefs.h:12
CAT_nKd_lin
#define CAT_nKd_lin
Definition: l2prod.h:372
i
int i
Definition: decode_rs.h:71
CAT_Kd_lee
#define CAT_Kd_lee
Definition: l2prod.h:108
Kd_lee
void Kd_lee(l2str *l2rec, int band, float *Kd)
spectral diffuse attenuation using Lee, et. (2005, 2016)
Definition: get_Kd.c:559
CAT_Kd_rhos
#define CAT_Kd_rhos
Definition: l2prod.h:322
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
a
PGE01 indicating that PGE02 PGE01 V6 for and PGE01 V2 for MOD03 were used to produce the granule By convention adopted in all MODIS Terra PGE02 code versions are The fourth digit of the PGE02 version denotes the LUT version used to produce the granule The source of the metadata environment variable ProcessingCenter was changed from a QA LUT value to the Process Configuration A sign used in error in the second order term was changed to a
Definition: HISTORY.txt:424
Kd490_unc
void Kd490_unc(l2str *l2rec, int32_t ip)
Definition: get_Kd.c:1323
npix
int npix
Definition: get_cmp.c:28
Zhl_morel
void Zhl_morel(l2str *l2rec, float *Zhl)
Definition: get_Kd.c:444
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
Kd490_KD2
void Kd490_KD2(l2str *l2rec, float *Kd)
Definition: get_Kd.c:29


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