OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
carder.c
Go to the documentation of this file.
1 /* =========================================================================
2  *
3  * Calculates aph675 and ag400 algebraically from Rrs model equations.
4  * chl is then calculated from aph675.
5  *
6  * Original Authors:
7  * Kendall Carder, F. R. Chen, and Steve Hawes
8  * University of South Florida
9  * Department of Marine Science
10  * 140 7th Avenue South
11  * St. Petersburg, FL 33701
12  * kcarder@monty.marine.usf.edu
13  *
14  * References:
15  * - Carder et al., Semi-Analytic MODIS Algorithms for Chlorophyll a and
16  * Absorption with Bio-Optical Domains Based on Nitrate-Depletion
17  * Temperatures, JGR MODIS special issue, 1998.
18  * - Carder et al., Reflectance model for quantifying chlorophyll _a_
19  * in the presence of productivity degradation products, JGR, 96(C11),
20  * 20599-20611, 1991.
21  * - Lee et al., Model for the interpretation of hyperspectral remote-
22  * sensing reflectance, Appl. Opt, 33(24), 5721, 1994.
23  *
24  * Modifications:
25  * 18 May 2004, B. Franz
26  * Generalized for use in MSl12. Pass in sensor wavelengths. Compute
27  * aw and bbw at input wavelengths. Add use of NDT for packaging switch
28  * based on modis_chl-1.2.c code. Code simplification. Enhanced flag
29  * tracking. General reorganization for efficiency.
30  *
31  * ========================================================================= */
32 
33 #include <math.h>
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include "l12_proto.h"
38 
39 #define N_PRM_HEADERS 35 /* leading comment lines in PARAM_FILE */
40 #define NX 32 /* elements minus one in aph675 array. */
41 #define N_ITER 5 /* iterations in bisection loop. 2^(N_ITER) = NX */
42 #define APH_MIN 0.0001 /* minimum aph675 value */
43 #define CARDER_FAIL -1 /* value of products returned if failure */
44 
45 /* flag bit settings */
46 #define CFLAG_FAIL 0x0001 /* failure (bad inputs) */
47 #define CFLAG_DEFAULT 0x0002 /* using default (empirical) algorithm*/
48 #define CFLAG_BLEND 0x0004 /* blending with default algorithm */
49 
50 #define CFLAG_UNPKG 0x0010 /* unpackaged case */
51 #define CFLAG_PKG 0x0020 /* packaged case */
52 #define CFLAG_HIPKG 0x0040 /* highly packaged case */
53 #define CFLAG_GLOBAL 0x0080 /* global case */
54 
55 #define CFLAG_LO412 0x0100 /* low input 412 reflectance */
56 #define CFLAG_LO555 0x0200 /* low input 55x reflectance */
57 #define CFLAG_CHLINC 0x0400 /* large inconsistency with default */
58 
59 static int pk; /* packaging parameter */
60 static int band1, band2, band3, band4; /* band indices to make Rrs ratios */
61 static int lam[6]; /* seawifs wavelengths */
62 static int bb_denom; /* toggle bb into denom of Rrs eqn */
63 static double bbw[6]; /* water backscattering coefs. */
64 static double aw[6]; /* water absorption coefs */
65 static double a0[6], a1[6], a2[6], a3[6]; /* aph model prms */
66 static double ga0[6], ga3[6]; /* aph model prms global */
67 static double upa0[6], upa3[6]; /* aph model prms unpackaged */
68 static double pa0[6], pa3[6]; /* aph model prms packaged */
69 static double hpa0[6], hpa1[6], hpa2[6], hpa3[6]; /* aph model prms hipackaged*/
70 static double s; /* ag exponential slope prm */
71 static double x0, x1; /* regr. coefs for x */
72 static double y_0, y_1; /* regr. coefs for y */
73 static double c0, c1, c2, c3; /* coefs for 490:555 chl alg */
74 static double gc0, gc1, gc2, gc3; /* coefs for 490:555 chl alg global */
75 static double upc0, upc1, upc2, upc3; /* coefs for 490:555 chl alg unpackaged*/
76 static double pc0, pc1, pc2, pc3; /* coefs for 490:555 chl alg packaged */
77 static double hpc0, hpc1, hpc2, hpc3; /* coefs for 490:555 chl alg hipackaged*/
78 static double p0, p1, p2; /* coefs for chl=fn(aph675) */
79 static double gp0, gp1, gp2; /* coefs for chl=fn(aph675) global */
80 static double upp0, upp1, upp2; /* coefs for chl=fn(aph675) unpackaged*/
81 static double pp0, pp1, pp2; /* coefs for chl=fn(aph675) packaged */
82 static double hpp0, hpp1, hpp2; /* coefs for chl=fn(aph675) hipackaged*/
83 static double d1, d2, d3; /* filter table parameters */
84 static double low_412_thresh, low_555_thresh, chl_inconsistent_thresh;
85 static double aph_hi, aph_lo; /* boundaries of aph675 array */
86 static double delta;
87 static double aphb1, aphb2, aphb3, aphb4, fmid;
88 static double tx[NX + 1];
89 static double arr0, arr1;
90 
91 /* function declarations */
92 
93 static void get_params(void);
94 static double get_aph(int band, double aph675);
95 
96 
97 /* ========================================================================= */
98 
99 /* ========================================================================= */
100 int carder_model(int32_t nbands, double *rrs, double *rrswl, float sst, float ndt,
101  double *atot, double *aph,
102  double *adg, double *bb, double *bbp,
103  double *aphmPtr, double *adgmPtr, double *chlPtr,
104  int16 *flags) {
105 
106  double br15, br25, br35, exponent, laph;
107  double adg_def, aph_def, chl_def;
108  double abr15, abr25, abr35;
109  double x, y, bb1, bb2, bb3, bb4; /* backscattering prms */
110  double r12, r34; /* Rrs/bb ratios */
111  double g12, g34, f0, f1, f2, f3, f4, wph;
112  double funk_lo, funk_hi, flo, fhi;
113  double wt;
114  double chl_mod, aph_mod, adg_mod; /* modeled chl,aph675,and ag400 */
115  double tchl_mod, taph_mod, tadg_mod; /* modeled chl,aph675,and ag400 (temp) */
116 
117  double weit; /* package switching */
118  double uptemp, gtemp, ptemp;
119  double sw0lo, sw1lo, sw2lo, sw3lo;
120  double dtsw1, dtsw2, dtsw3;
121  int it, pktran[2];
122 
123  int i, xlo, xmid, xhi;
124  static double cpa1[6], cpa2[6]; /* copies of a1 and a2 */
125  static int firstCall = 1;
126 
127  if (firstCall) {
128 
129  firstCall = 0;
130 
131  /* read model parameters */
132  get_params();
133 
134  /* copy the a1 and a2 arrays */
135  for (i = 0; i < 6; i++) cpa1[i] = a1[i];
136  for (i = 0; i < 6; i++) cpa2[i] = a2[i];
137 
138  /* create aph675 array, 0.0001-aph_hi (1/m) */
139  arr0 = log10(aph_lo);
140  arr1 = log10(aph_hi) - arr0;
141  for (i = 0; i <= NX; i++)
142  tx[i] = pow(10.0, arr0 + arr1 * (double) i / (double) NX);
143 
144  /* over-ride parameter inputs with sensor-band-specific aw, bbw */
145  /* using 10-nm square bandpass average to be consistent with Rrs */
146  for (i = 0; i < 6; i++) {
147  if (rrswl[i] > 0.0) {
148  aw [i] = (double) aw_spectra(rrswl[i], BANDW);
149  bbw[i] = (double) bbw_spectra(rrswl[i], BANDW);
150  } else {
151  aw [i] = (double) 1.0;
152  bbw[i] = (double) 1.0;
153  }
154  }
155  }
156 
157  /* initialize */
158  for (i = 0; i < nbands; i++) {
159  bb [i] = CARDER_FAIL;
160  atot[i] = CARDER_FAIL;
161  aph [i] = CARDER_FAIL;
162  adg [i] = CARDER_FAIL;
163  }
164  *chlPtr = CARDER_FAIL;
165  *adgmPtr = CARDER_FAIL;
166  *aphmPtr = CARDER_FAIL;
167  *flags = 0x0000;
168 
169  /* if we have any negative reflectances, just bail */
170  if (rrs[0] <= 0. || rrs[1] <= 0. || rrs[2] <= 0. || rrs[4] <= 0.) {
171  *flags |= CFLAG_FAIL;
172  return 1;
173  }
174 
175  /* set some quality flags */
176  if (rrs[0] < low_412_thresh) {
177  *flags |= CFLAG_LO412;
178  }
179  if (rrs[4] < low_555_thresh) {
180  *flags |= CFLAG_LO555;
181  }
182 
183  /* set up the limits based on NDT */
184  uptemp = ndt + 2.;
185  gtemp = ndt + 0.8;
186  ptemp = ndt - 1.;
187 
188  sw0lo = uptemp + 1.;
189  sw1lo = gtemp + (uptemp - gtemp) / 2.;
190  sw2lo = ptemp + (gtemp - ptemp) / 2.;
191  sw3lo = ptemp - 1.;
192 
193  dtsw1 = sw0lo - sw1lo;
194  dtsw2 = sw1lo - sw2lo;
195  dtsw3 = sw2lo - sw3lo;
196 
197  /* sw = 0 */
198  if (sst >= sw0lo) {
199  weit = 1.;
200  pktran[0] = 1;
201  pktran[1] = 1;
202  }
203 
204  /* sw = 1 */
205  if (sst >= sw1lo && sst < sw0lo) {
206  weit = (sst - sw1lo) / dtsw1;
207  pktran[0] = 0;
208  pktran[1] = 1;
209  }
210 
211  /* sw = 2 */
212  if (sst >= sw2lo && sst < sw1lo) {
213  weit = (sst - sw2lo) / dtsw2;
214  pktran[0] = 2;
215  pktran[1] = 0;
216  }
217 
218  /* sw = 3 */
219  if (sst >= sw3lo && sst < sw2lo) {
220  weit = (sst - sw3lo) / dtsw3;
221  pktran[0] = 3;
222  pktran[1] = 2;
223  }
224 
225  /* sw = 4 */
226  if (sst < sw3lo) {
227  weit = 1.;
228  pktran[0] = 3;
229  pktran[1] = 3;
230  }
231 
232 
233  /* calculate some bandratios */
234  br15 = rrs[0] / rrs[4];
235  br25 = rrs[1] / rrs[4];
236  br35 = rrs[2] / rrs[4];
237  abr15 = log10(br15);
238  abr25 = log10(br25);
239  abr35 = log10(br35);
240 
241  /* calculate default ag400 and aph675 */
242  exponent = MAX(MIN(-1.147 - 1.963 * abr15 - 1.01 * abr15 * abr15 + 0.856 * abr25 + 1.702 * abr25*abr25, 10.0), -10.0);
243  adg_def = 1.5 * pow(10.0, exponent);
244  exponent = MAX(MIN(-0.919 + 1.037 * abr25 - 0.407 * abr25 * abr25 - 3.531 * abr35 + 1.579 * abr35*abr35, 10.0), -10.0);
245  aph_def = MAX((pow(10.0, exponent) - 0.008) / 3.05, APH_MIN);
246 
247  /* calculate bb */
248  if (rrs[4] > 0.001)
249  x = MAX(x0 + x1 * rrs[4], 0.0);
250  else
251  x = MAX(x0 + x1 * 0.001, 0.0);
252 
253  y = MAX(y_0 + y_1 * rrs[1] / rrs[2], 0.0);
254 
255  for (i = 0; i < 6; i++) {
256  bbp[i] = x * pow(551. / rrswl[i], y);
257  bb [i] = bbw[i] + bbp[i];
258  }
259 
260  /* calc. f0, f1, f2, f3, f4; coefs. for funk() */
261  bb1 = bb[band1];
262  bb2 = bb[band2];
263  bb3 = bb[band3];
264  bb4 = bb[band4];
265  r12 = (rrs[band1] / bb1) / (rrs[band2] / bb2);
266  r34 = (rrs[band3] / bb3) / (rrs[band4] / bb4);
267  g12 = r12 * exp(-s * (rrswl[band1] - 400.)) - exp(-s * (rrswl[band2] - 400.));
268  g34 = r34 * exp(-s * (rrswl[band3] - 400.)) - exp(-s * (rrswl[band4] - 400.));
269  if (bb_denom) {
270  f0 = g12 * (aw[band4] + bb4 - r34 * (aw[band3] + bb3))
271  - g34 * (aw[band2] + bb2 - r12 * (aw[band1] + bb1));
272  } else {
273  f0 = g12 * (aw[band4] - r34 * aw[band3])
274  - g34 * (aw[band2] - r12 * aw[band1]);
275  }
276  f1 = g34*r12;
277  f2 = -g34;
278  f3 = -g12*r34;
279  f4 = g12;
280 
281 
282  /* Loop through twice, to blend packaging effects */
283  for (it = 0; it <= 1; it++) {
284 
285  pk = pktran[it];
286 
287  switch (pk) {
288  case 1: /* unpackaged */
289  a0[0] = upa0[0];
290  a0[1] = upa0[1];
291  a0[2] = upa0[2];
292  a0[3] = upa0[3];
293  a0[4] = upa0[4];
294  a0[5] = upa0[5];
295  a1[0] = cpa1[0];
296  a1[1] = cpa1[1];
297  a1[2] = cpa1[2];
298  a1[3] = cpa1[3];
299  a1[4] = cpa1[4];
300  a1[5] = cpa1[5];
301  a2[0] = cpa2[0];
302  a2[1] = cpa2[1];
303  a2[2] = cpa2[2];
304  a2[3] = cpa2[3];
305  a2[4] = cpa2[4];
306  a2[5] = cpa2[5];
307  a3[0] = upa3[0];
308  a3[1] = upa3[1];
309  a3[2] = upa3[2];
310  a3[3] = upa3[3];
311  a3[4] = upa3[4];
312  a3[5] = upa3[5];
313  c0 = upc0;
314  c1 = upc1;
315  c2 = upc2;
316  c3 = upc3;
317  p0 = upp0;
318  p1 = upp1;
319  p2 = upp2;
320  *flags |= CFLAG_UNPKG;
321  break;
322  case 2: /* packaged */
323  a0[0] = pa0[0];
324  a0[1] = pa0[1];
325  a0[2] = pa0[2];
326  a0[3] = pa0[3];
327  a0[4] = pa0[4];
328  a0[5] = pa0[5];
329  a1[0] = cpa1[0];
330  a1[1] = cpa1[1];
331  a1[2] = cpa1[2];
332  a1[3] = cpa1[3];
333  a1[4] = cpa1[4];
334  a1[5] = cpa1[5];
335  a2[0] = cpa2[0];
336  a2[1] = cpa2[1];
337  a2[2] = cpa2[2];
338  a2[3] = cpa2[3];
339  a2[4] = cpa2[4];
340  a2[5] = cpa2[5];
341  a3[0] = pa3[0];
342  a3[1] = pa3[1];
343  a3[2] = pa3[2];
344  a3[3] = pa3[3];
345  a3[4] = pa3[4];
346  a3[5] = pa3[5];
347  c0 = pc0;
348  c1 = pc1;
349  c2 = pc2;
350  c3 = pc3;
351  p0 = pp0;
352  p1 = pp1;
353  p2 = pp2;
354  *flags |= CFLAG_PKG;
355  break;
356  case 3: /*hipackaged */
357  a0[0] = hpa0[0];
358  a0[1] = hpa0[1];
359  a0[2] = hpa0[2];
360  a0[3] = hpa0[3];
361  a0[4] = hpa0[4];
362  a0[5] = hpa0[5];
363  a1[0] = hpa1[0];
364  a1[1] = hpa1[1];
365  a1[2] = hpa1[2];
366  a1[3] = hpa1[3];
367  a1[4] = hpa1[4];
368  a1[5] = hpa1[5];
369  a2[0] = hpa2[0];
370  a2[1] = hpa2[1];
371  a2[2] = hpa2[2];
372  a2[3] = hpa2[3];
373  a2[4] = hpa2[4];
374  a2[5] = hpa2[5];
375  a3[0] = hpa3[0];
376  a3[1] = hpa3[1];
377  a3[2] = hpa3[2];
378  a3[3] = hpa3[3];
379  a3[4] = hpa3[4];
380  a3[5] = hpa3[5];
381  c0 = hpc0;
382  c1 = hpc1;
383  c2 = hpc2;
384  c3 = hpc3;
385  p0 = hpp0;
386  p1 = hpp1;
387  p2 = hpp2;
388  *flags |= CFLAG_HIPKG;
389  break;
390  default: /* global */
391  a0[0] = ga0[0];
392  a0[1] = ga0[1];
393  a0[2] = ga0[2];
394  a0[3] = ga0[3];
395  a0[4] = ga0[4];
396  a0[5] = ga0[5];
397  a1[0] = cpa1[0];
398  a1[1] = cpa1[1];
399  a1[2] = cpa1[2];
400  a1[3] = cpa1[3];
401  a1[4] = cpa1[4];
402  a1[5] = cpa1[5];
403  a2[0] = cpa2[0];
404  a2[1] = cpa2[1];
405  a2[2] = cpa2[2];
406  a2[3] = cpa2[3];
407  a2[4] = cpa2[4];
408  a2[5] = cpa2[5];
409  a3[0] = ga3[0];
410  a3[1] = ga3[1];
411  a3[2] = ga3[2];
412  a3[3] = ga3[3];
413  a3[4] = ga3[4];
414  a3[5] = ga3[5];
415  c0 = gc0;
416  c1 = gc1;
417  c2 = gc2;
418  c3 = gc3;
419  p0 = gp0;
420  p1 = gp1;
421  p2 = gp2;
422  *flags |= CFLAG_GLOBAL;
423  break;
424  }
425 
426 
427  /* calculate default chl, note coeffs depend on packaging */
428  exponent = MAX(MIN(c0 + c1 * abr35 + c2 * abr35 * abr35 + c3 * abr35 * abr35*abr35, 10.0), -10.0);
429  chl_def = pow(10.0, exponent);
430 
431 
432  /* solve model */
433  /* ----------- */
434 
435  aphb1 = get_aph(band1, tx[NX]);
436  aphb2 = get_aph(band2, tx[NX]);
437  aphb3 = get_aph(band3, tx[NX]);
438  aphb4 = get_aph(band4, tx[NX]);
439  funk_hi = f0 + f1 * aphb1 + f2 * aphb2 + f3 * aphb3 + f4*aphb4;
440 
441  aphb1 = get_aph(band1, tx[0]);
442  aphb2 = get_aph(band2, tx[0]);
443  aphb3 = get_aph(band3, tx[0]);
444  aphb4 = get_aph(band4, tx[0]);
445  funk_lo = f0 + f1 * aphb1 + f2 * aphb2 + f3 * aphb3 + f4*aphb4;
446 
447  if (funk_lo > 0.0 && funk_hi < 0.0) {
448 
449  /* solve for aph675 via bisection */
450 
451  xlo = 0;
452  xhi = NX;
453 
454  for (i = 1; i <= N_ITER; i++) {
455 
456  xmid = (xlo + xhi) / 2;
457  aphb1 = get_aph(band1, tx[xmid]);
458  aphb2 = get_aph(band2, tx[xmid]);
459  aphb3 = get_aph(band3, tx[xmid]);
460  aphb4 = get_aph(band4, tx[xmid]);
461  fmid = f0 + f1 * aphb1 + f2 * aphb2 + f3 * aphb3 + f4*aphb4;
462  if (fmid < 0.0)
463  xhi = xmid;
464  else
465  xlo = xmid;
466  }
467 
468  aphb1 = get_aph(band1, tx[xlo]);
469  aphb2 = get_aph(band2, tx[xlo]);
470  aphb3 = get_aph(band3, tx[xlo]);
471  aphb4 = get_aph(band4, tx[xlo]);
472  flo = f0 + f1 * aphb1 + f2 * aphb2 + f3 * aphb3 + f4*aphb4;
473 
474  aphb1 = get_aph(band1, tx[xhi]);
475  aphb2 = get_aph(band2, tx[xhi]);
476  aphb3 = get_aph(band3, tx[xhi]);
477  aphb4 = get_aph(band4, tx[xhi]);
478  fhi = f0 + f1 * aphb1 + f2 * aphb2 + f3 * aphb3 + f4*aphb4;
479 
480  /* aph675; linear interp across funk to find zero */
481  taph_mod = tx[xlo] + (tx[xhi] - tx[xlo]) * flo / (flo - fhi);
482 
483  /* ag400 = fn(aph675) */
484  if (bb_denom) {
485  wph = aw[band4] + get_aph(band4, taph_mod) + bb4
486  - r34 * (aw[band3] + get_aph(band3, taph_mod) + bb3);
487  } else {
488  wph = aw[band4] + get_aph(band4, taph_mod)
489  - r34 * (aw[band3] + get_aph(band3, taph_mod));
490  }
491  tadg_mod = wph / g34;
492 
493  /* chlorophyll */
494  laph = log10(taph_mod);
495  exponent = p0 + p1 * laph + p2 * laph * laph;
496  tchl_mod = pow(10.0, exponent);
497 
498  if (fabs(log10(tchl_mod / chl_def)) > chl_inconsistent_thresh) {
499  *flags |= CFLAG_CHLINC;
500  }
501 
502  /* if aph_mod > aph_hi/2, blend with default alg. */
503  if (taph_mod > aph_hi / 2.) {
504  wt = (tx[NX] - taph_mod) / (tx[NX] - aph_hi / 2.);
505  tchl_mod = wt * tchl_mod + (1.0 - wt) * chl_def;
506  tadg_mod = wt * tadg_mod + (1.0 - wt) * adg_def;
507  taph_mod = wt * taph_mod + (1.0 - wt) * aph_def;
508  *flags |= CFLAG_BLEND;
509  }
510 
511  } else {
512 
513  /* if funk doesn't bracket zero, return defaults */
514 
515  tchl_mod = chl_def;
516  tadg_mod = adg_def;
517  taph_mod = aph_def;
518  *flags |= CFLAG_DEFAULT;
519  }
520 
521  /* weight results */
522  if (it == 0) {
523  chl_mod = tchl_mod;
524  adg_mod = tadg_mod;
525  aph_mod = taph_mod;
526  for (i = 0; i < 6; i++)
527  aph[i] = get_aph(i, aph_mod);
528  } else {
529  chl_mod = chl_mod * (1. - weit) + tchl_mod * weit;
530  adg_mod = adg_mod * (1. - weit) + tadg_mod * weit;
531  aph_mod = aph_mod * (1. - weit) + taph_mod * weit;
532  for (i = 0; i < 6; i++)
533  aph[i] = aph[i]*(1. - weit) + get_aph(i, aph_mod) * weit;
534  }
535  }
536 
537  /* Store results for output */
538  for (i = 0; i < 6; i++) {
539  adg [i] = adg_mod * exp(-s * (rrswl[i] - 400.0));
540  atot[i] = aw[i] + aph[i] + adg[i];
541  }
542 
543  *aphmPtr = aph_mod; /* model value aph(675) */
544  *adgmPtr = adg_mod; /* model value adg(400) */
545  *chlPtr = chl_mod;
546 
547  return 0;
548 
549 }
550 
551 /******************** FUNCTION get_params ******************************
552  *
553  * Reads in the the seawifs wavebands and Rrs model parameters to use
554  * in the algorithm.
555  *
556  *************************************************************************/
557 
558 static void get_params(void) {
559  int i;
560  char dummy[100], full_file[2048];
561  FILE *fin;
562  char *tmp_str;
563 
564  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
565  printf("OCDATAROOT environment variable is not defined.\n");
566  exit(1);
567  }
568 
569  strcpy(full_file, tmp_str);
570  strcat(full_file, "/common/carder.par");
571  fin = fopen(full_file, "r");
572  if (!fin) {
573  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, full_file);
574  exit(1);
575  }
576 
577  printf("\nLoading Carder parameters from file %s\n", full_file);
578 
579  for (i = 1; i <= N_PRM_HEADERS; i++)
580  fgets(dummy, 100 - 1, fin);
581  fscanf(fin, "%d\n", &pk);
582  fscanf(fin, "%lf %lf %lf %lf\n", &d1, &d2, &d3, &delta); /* filter table param*/
583  fscanf(fin, "%d %d %d %d\n", &band1, &band2, &band3, &band4);
584  for (i = 0; i < 6; i++)
585  fscanf(fin, "%d %lf %lf %lf %lf %lf %lf\n",
586  &lam[i], &bbw[i], &aw[i], &ga0[i], &a1[i], &a2[i], &ga3[i]);
587  fscanf(fin, "%lf\n", &s); /* ag exp. slope */
588  fscanf(fin, "%lf %lf\n", &x0, &x1); /* regr. coefs for x */
589  fscanf(fin, "%lf %lf\n", &y_0, &y_1); /* regr. coefs for y */
590  fscanf(fin, "%lf %lf %lf %lf\n", &gc0, &gc1, &gc2, &gc3); /* coefs for r35 chl*/
591  fscanf(fin, "%lf %lf %lf\n", &gp0, &gp1, &gp2); /* coefs for chl-aph675 */
592  fscanf(fin, "%lf %lf %lf\n", &low_412_thresh, &low_555_thresh,
593  &chl_inconsistent_thresh);
594  fscanf(fin, "%lf %lf\n", &aph_lo, &aph_hi);
595  fscanf(fin, "%d\n", &bb_denom);
596  fscanf(fin, "%lf %lf %lf %lf %lf %lf\n", &upa0[0], &upa0[1],
597  &upa0[2], &upa0[3], &upa0[4], &upa0[5]);
598  fscanf(fin, "%lf %lf %lf %lf %lf %lf\n", &upa3[0], &upa3[1],
599  &upa3[2], &upa3[3], &upa3[4], &upa3[5]);
600  fscanf(fin, "%lf %lf %lf %lf %lf %lf\n", &pa0[0], &pa0[1],
601  &pa0[2], &pa0[3], &pa0[4], &pa0[5]);
602  fscanf(fin, "%lf %lf %lf %lf %lf %lf\n", &pa3[0], &pa3[1],
603  &pa3[2], &pa3[3], &pa3[4], &pa3[5]);
604  fscanf(fin, "%lf %lf %lf %lf %lf %lf\n", &hpa0[0], &hpa0[1],
605  &hpa0[2], &hpa0[3], &hpa0[4], &hpa0[5]);
606  fscanf(fin, "%lf %lf %lf %lf %lf %lf\n", &hpa1[0], &hpa1[1],
607  &hpa1[2], &hpa1[3], &hpa1[4], &hpa1[5]);
608  fscanf(fin, "%lf %lf %lf %lf %lf %lf\n", &hpa2[0], &hpa2[1],
609  &hpa2[2], &hpa2[3], &hpa2[4], &hpa2[5]);
610  fscanf(fin, "%lf %lf %lf %lf %lf %lf\n", &hpa3[0], &hpa3[1],
611  &hpa3[2], &hpa3[3], &hpa3[4], &hpa3[5]);
612  fscanf(fin, "%lf %lf %lf %lf\n", &upc0, &upc1, &upc2, &upc3);
613  fscanf(fin, "%lf %lf %lf %lf\n", &pc0, &pc1, &pc2, &pc3);
614  fscanf(fin, "%lf %lf %lf %lf\n", &hpc0, &hpc1, &hpc2, &hpc3);
615  fscanf(fin, "%lf %lf %lf\n", &upp0, &upp1, &upp2);
616  fscanf(fin, "%lf %lf %lf\n", &pp0, &pp1, &pp2);
617  fscanf(fin, "%lf %lf %lf\n", &hpp0, &hpp1, &hpp2);
618  fclose(fin);
619 }
620 
621 /******************** FUNCTION aph *************************************
622  *
623  * Returns a_phi at a given waveband as a function of aph675.
624  *
625  *************************************************************************/
626 
627 static double get_aph(int band, double aph675) {
628  if (aph675 == -1.0) {
629  return (-1.0);
630  } else {
631  return a0[band] * exp(a1[band] * tanh(a2[band] * log(aph675 / a3[band]))) * aph675;
632  }
633 }
634 
635 
636 
637 
638 /* ----------------------------------------------------------------------- */
639 /* get_ndt() - read and interpolate NDT map */
640 
641 /* ----------------------------------------------------------------------- */
642 static float get_ndt(float lon, float lat) {
643 
644 #define NDTNX 2048
645 #define NDTNY 1024
646 
647  static int firstCall = 1;
648 
649  typedef float map_t[NDTNX];
650  static map_t* map;
651  static float slope;
652  static float offset;
653 
654  int i, j;
655 
656  if (firstCall) {
657 
658  char *tmp_str;
659  int32 sd_id;
660  int32 sds_id;
661  int32 status;
662  int32 sds_index;
663  int32 rank;
664  int32 nt;
665  int32 dims[H4_MAX_VAR_DIMS];
666  int32 nattrs;
667  int32 start[2];
668  int32 edges[2];
669  char name[H4_MAX_NC_NAME];
670  char sdsname[] = "ndt";
671  char file[FILENAME_MAX];
672 
673  firstCall = 0;
674 
675  // allocate map array
676  map = (map_t*) allocateMemory(NDTNY * sizeof (map_t), "ndt map");
677 
678 
679  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
680  printf("OCDATAROOT environment variable is not defined.\n");
681  exit(1);
682  }
683  strcpy(file, tmp_str);
684  strcat(file, "/common/ndt.hdf");
685 
686  printf("\nLoading nitrate depletion temperature map from %s\n", file);
687 
688  /* Open the file and initiate the SD interface */
689  sd_id = SDstart(file, DFACC_RDONLY);
690  if (sd_id == -1) {
691  printf("-E- %s: Error opening file %s.\n",
692  __FILE__, file);
693  exit(1);
694  }
695 
696  /* Get the SDS index */
697  sds_index = SDnametoindex(sd_id, sdsname);
698  if (sds_index == -1) {
699  printf("-E- %s: Error seeking %s SDS from %s.\n",
700  __FILE__, sdsname, file);
701  exit(1);
702  }
703 
704  /* Select the SDS */
705  sds_id = SDselect(sd_id, sds_index);
706 
707  /* Verify the characteristics of the array */
708  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
709  if (status != 0) {
710  printf("-E- %s: Error reading SDS info for %s from %s.\n",
711  __FILE__, sdsname, file);
712  exit(1);
713  }
714  if (dims[0] != NDTNY || dims[1] != NDTNX) {
715  printf("-E- %s: Dimension mis-match on %s array from %s.\n",
716  __FILE__, sdsname, file);
717  printf(" Expecting %d x %d\n", NDTNX, NDTNY);
718  printf(" Reading %d x %d\n", dims[1], dims[0]);
719  exit(1);
720  }
721 
722  start[0] = 0; /* row offset */
723  start[1] = 0; /* col offset */
724  edges[0] = NDTNY; /* row count */
725  edges[1] = NDTNX; /* col count */
726 
727  status = SDreaddata(sds_id, start, NULL, edges, (VOIDP) map);
728  if (status != 0) {
729  printf("-E- %s: Error reading SDS %s from %s.\n",
730  __FILE__, sdsname, file);
731  exit(1);
732  }
733 
734  /* Read the scaling info */
735  status = SDreadattr(sds_id, SDfindattr(sds_id, "slope"), (VOIDP) & slope);
736  if (status != 0) {
737  printf("-E- %s: Error reading slope attribute for SDS %s from %s.\n",
738  __FILE__, sdsname, file);
739  exit(1);
740  }
741  status = SDreadattr(sds_id, SDfindattr(sds_id, "intercept"), (VOIDP) & offset);
742  if (status != 0) {
743  printf("-E- %s: Error reading intercept attribute for SDS %s from %s.\n",
744  __FILE__, sdsname, file);
745  exit(1);
746  }
747 
748  /* Terminate access to the array */
749  status = SDendaccess(sds_id);
750 
751  /* Terminate access to the SD interface and close the file */
752  status = SDend(sd_id);
753  }
754 
755  /* No reason for perfection here */
756 
757  i = MIN(MAX((int) ((lon + 180.) / 360.0 * NDTNX), 0), NDTNX - 1);
758  j = MIN(MAX((int) ((lat + 90.) / 180.0 * NDTNY), 0), NDTNY - 1);
759 
760  return (map[j][i] * slope + offset);
761 }
762 
763 
764 /* ----------------------------------------------------------------------- */
765 /* MSl12 interface to Carder model (from carder_iop.c) */
766 /* ----------------------------------------------------------------------- */
767 
768 #include <stdlib.h>
769 #include <math.h>
770 #include "l12_proto.h"
771 
772 /* These hold 2-D arrays of BIP for each of the products calculated by carder. */
773 
774 static double *atot;
775 static double *aph, *aphm;
776 static double *adg, *adgm;
777 static double *bbp, *bb;
778 static double *chl;
779 static int16 *flagm;
780 static int CarderRecNum = -99;
781 
782 /* Allocate private arrays for a single scan line */
783 
784 void alloc_carder(int32_t npix, int32_t nbands) {
785  atot = (double*) malloc(npix * nbands * sizeof (double));
786  aph = (double*) malloc(npix * nbands * sizeof (double));
787  adg = (double*) malloc(npix * nbands * sizeof (double));
788  bb = (double*) malloc(npix * nbands * sizeof (double));
789  bbp = (double*) malloc(npix * nbands * sizeof (double));
790  chl = (double*) malloc(npix * sizeof (double));
791  aphm = (double*) malloc(npix * sizeof (double));
792  adgm = (double*) malloc(npix * sizeof (double));
793  flagm = (int16*) malloc(npix * sizeof (int16));
794 }
795 
796 /* Determine if we have already run carder's algorithm for this scan line */
797 
798 int carder_ran(int recnum) {
799 
800  if (recnum == CarderRecNum)
801  return 1;
802  else
803  return 0;
804 
805 }
806 
807 /* Set-up and run the Carder model for each pixel in this scan */
808 
809 void carder(l2str *l2rec) {
810  static double *rrswl;
811  static int *indx;
812  static int firstCall = 1;
813 
814  double *rrs;
815  double *aph1;
816  double *adg1;
817  double *atot1;
818 
819  double *bb1;
820  double *bbp1;
821  double aph_mod; /* aph(675) */
822  double adg_mod; /* adg(400) */
823  double chl_mod;
824  int16 flags;
825  int ip, ic;
826  int ib;
827  int32_t ipb;
828  float ndt;
829  float sst;
830  int status, nbands;
831 
832  l1str *l1rec = l2rec->l1rec;
833  nbands = l1rec->l1file->nbands;
834 
835 
836  if (firstCall) {
837  firstCall = 0;
838  if ((rrswl = (double *) calloc(nbands, sizeof (double))) == NULL) {
839  printf("-E- : Error allocating memory to rrswl\n");
840  exit(FATAL_ERROR);
841  }
842  if ((indx = (int *) calloc(nbands, sizeof (int))) == NULL) {
843  printf("-E- : Error allocating memory to indx\n");
844  exit(FATAL_ERROR);
845  }
846 
847  alloc_carder(l1rec->npix, nbands);
848  indx[0] = bindex_get(412);
849  indx[1] = bindex_get(443);
850  indx[2] = bindex_get(490);
851  indx[3] = bindex_get(510);
852  if (indx[3] < 0) indx[3] = bindex_get(520);
853  if (indx[3] < 0) indx[3] = bindex_get(530);
854  indx[4] = bindex_get(551);
855  indx[5] = bindex_get(670);
856  for (ib = 0; ib < 6; ib++) {
857  rrswl[ib] = (double) l1rec->l1file->fwave[indx[ib]];
858  printf("Carder wavelength %d is %f\n", ib, rrswl[ib]);
859  }
860  }
861 
862  if ((rrs = (double *) calloc(nbands, sizeof (double))) == NULL) {
863  printf("-E- : Error allocating memory to rrs\n");
864  exit(FATAL_ERROR);
865  }
866  if ((aph1 = (double *) calloc(nbands, sizeof (double))) == NULL) {
867  printf("-E- : Error allocating memory to aph1\n");
868  exit(FATAL_ERROR);
869  }
870  if ((adg1 = (double *) calloc(nbands, sizeof (double))) == NULL) {
871  printf("-E- : Error allocating memory to adg1\n");
872  exit(FATAL_ERROR);
873  }
874  if ((atot1 = (double *) calloc(nbands, sizeof (double))) == NULL) {
875  printf("-E- : Error allocating memory to atot1\n");
876  exit(FATAL_ERROR);
877  }
878  if ((bb1 = (double *) calloc(nbands, sizeof (double))) == NULL) {
879  printf("-E- : Error allocating memory to bb1\n");
880  exit(FATAL_ERROR);
881  }
882  if ((bbp1 = (double *) calloc(nbands, sizeof (double))) == NULL) {
883  printf("-E- : Error allocating memory to bbp1\n");
884  exit(FATAL_ERROR);
885  }
886 
887  for (ip = 0; ip < l1rec->npix; ip++) {
888 
889  ipb = ip*nbands;
890 
891  /* grab the remote-sensing reflectances for this pixel */
892  for (ic = 0; ic < 6; ic++)
893  rrs[ic] = (double) (l2rec->Rrs[ipb + indx[ic]]);
894 
895  /* Use retrieved sst, OR reference sst if no retrieval */
896  if (l2rec->sst != NULL && l2rec->sst[ip] > -2.0)
897  sst = l2rec->sst[ip];
898  else
899  sst = l1rec->sstref[ip];
900 
901  /* grab ndt for this lon/lat */
902  ndt = get_ndt(l1rec->lon[ip], l1rec->lat[ip]);
903 
904  /* run the model */
905  status = carder_model(nbands, rrs, rrswl, sst, ndt, atot1, aph1, adg1, bb1, bbp1, &aph_mod, &adg_mod, &chl_mod, &flags);
906 
907  /* store results of this pixel into static global arrays */
908  if (status == 0) {
909  for (ic = 0; ic < 6; ic++) {
910  atot[ipb + indx[ic]] = atot1[ic];
911  aph [ipb + indx[ic]] = aph1 [ic];
912  adg [ipb + indx[ic]] = adg1 [ic];
913  bb [ipb + indx[ic]] = bb1 [ic];
914  bbp [ipb + indx[ic]] = bbp1 [ic];
915  }
916  flagm[ip] = flags;
917  adgm [ip] = adg_mod;
918  aphm [ip] = aph_mod;
919  chl [ip] = chl_mod;
920  } else {
921  for (ic = 0; ic < 6; ic++) {
922  atot[ipb + indx[ic]] = CARDER_FAIL;
923  aph [ipb + indx[ic]] = CARDER_FAIL;
924  adg [ipb + indx[ic]] = CARDER_FAIL;
925  bb [ipb + indx[ic]] = CARDER_FAIL;
926  bbp [ipb + indx[ic]] = CARDER_FAIL;
927  }
928  flagm[ip] = flags;
929  adgm [ip] = CARDER_FAIL;
930  aphm [ip] = CARDER_FAIL;
931  chl [ip] = CARDER_FAIL;
932  l1rec->flags[ip] |= PRODFAIL;
933  }
934  }
935 
936  CarderRecNum = l1rec->iscan;
937 
938  free(rrs);
939  free(aph1);
940  free(adg1);
941  free(atot1);
942  free(bb1);
943  free(bbp1);
944 }
945 
946 /* Interface to l2_hdf_generic() to return Carder flags */
947 
948 int16 *get_flags_carder(l2str *l2rec) {
949  if (!carder_ran(l2rec->l1rec->iscan))
950  carder(l2rec);
951 
952  return (flagm);
953 }
954 
955 /* Interface to l2_hdf_generic() to return Carder products */
956 
957 void get_carder(l2str *l2rec, l2prodstr *p, float prod[]) {
958  int band = p->prod_ix;
959  int prodID = p->cat_ix;
960  int ip, ipb;
961 
962  if (!carder_ran(l2rec->l1rec->iscan))
963  carder(l2rec);
964 
965  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
966 
967  ipb = ip * l2rec->l1rec->l1file->nbands + band;
968 
969  switch (prodID) {
970 
971  case CAT_chl_carder:
972  prod[ip] = (float) chl[ip];
973  break;
974 
975  case CAT_adg_carder:
976  if (band < 0)
977  prod[ip] = (float) adgm[ip];
978  else
979  prod[ip] = (float) adg[ipb];
980  break;
981 
982  case CAT_aph_carder:
983  if (band < 0)
984  prod[ip] = (float) aphm[ip];
985  else
986  prod[ip] = (float) aph[ipb];
987  break;
988 
989  case CAT_a_carder:
990  prod[ip] = (float) atot[ipb];
991  break;
992 
993  case CAT_bbp_carder:
994  prod[ip] = (float) bbp[ipb];
995  break;
996 
997  case CAT_bb_carder:
998  prod[ip] = (float) bb[ipb];
999  break;
1000 
1001  default:
1002  printf("-E- %s line %d : erroneous product ID %d passed to get_carder().\n",
1003  __FILE__, __LINE__, prodID);
1004  exit(1);
1005  }
1006  }
1007 
1008  return;
1009 }
1010 
1011 /* Interface to convl12() to return Carder iops */
1012 
1013 void iops_carder(l2str *l2rec) {
1014  int32_t ib, ip, ipb;
1015  int32_t nbands = l2rec->l1rec->l1file->nbands;
1016 
1017  if (!carder_ran(l2rec->l1rec->iscan))
1018  carder(l2rec);
1019 
1020  for (ip = 0; ip < l2rec->l1rec->npix; ip++) for (ib = 0; ib < nbands; ib++) {
1021  ipb = ip * nbands + ib;
1022  l2rec->a [ipb] = (float) atot[ipb];
1023  l2rec->bb[ipb] = (float) bb [ipb];
1024  }
1025 
1026  return;
1027 }
1028 
1029 
1030 
1031 /* ========================================================================= */
1032 
1033 /* ========================================================================= */
1034 int carder_empirical(float *rrs, float sst, float ndt, float *chlPtr) {
1035  static int firstCall = 1;
1036 
1037  float exponent;
1038  float abr35;
1039  float tchl;
1040  float chl;
1041  float weit; /* package switching */
1042 
1043  float uptemp, gtemp, ptemp;
1044  float sw0lo, sw1lo, sw2lo, sw3lo;
1045  float dtsw1, dtsw2, dtsw3;
1046  int it, pktran[2];
1047 
1048  if (firstCall) {
1049 
1050  firstCall = 0;
1051 
1052  /* read model parameters */
1053  get_params();
1054  }
1055 
1056  /* initialize */
1057  *chlPtr = CARDER_FAIL;
1058 
1059  /* if we have any negative reflectances, just bail */
1060  if (rrs[2] <= 0. || rrs[4] <= 0.) {
1061  return 1;
1062  }
1063 
1064  /* calculate bandratio */
1065  abr35 = log10(rrs[2] / rrs[4]);
1066 
1067  /* set up the limits based on NDT */
1068  uptemp = ndt + 2.;
1069  gtemp = ndt + 0.8;
1070  ptemp = ndt - 1.;
1071 
1072  sw0lo = uptemp + 1.;
1073  sw1lo = gtemp + (uptemp - gtemp) / 2.;
1074  sw2lo = ptemp + (gtemp - ptemp) / 2.;
1075  sw3lo = ptemp - 1.;
1076 
1077  dtsw1 = sw0lo - sw1lo;
1078  dtsw2 = sw1lo - sw2lo;
1079  dtsw3 = sw2lo - sw3lo;
1080 
1081  /* sw = 0 */
1082  if (sst >= sw0lo) {
1083  weit = 1.;
1084  pktran[0] = 1;
1085  pktran[1] = 1;
1086  }
1087 
1088  /* sw = 1 */
1089  if (sst >= sw1lo && sst < sw0lo) {
1090  weit = (sst - sw1lo) / dtsw1;
1091  pktran[0] = 0;
1092  pktran[1] = 1;
1093  }
1094 
1095  /* sw = 2 */
1096  if (sst >= sw2lo && sst < sw1lo) {
1097  weit = (sst - sw2lo) / dtsw2;
1098  pktran[0] = 2;
1099  pktran[1] = 0;
1100  }
1101 
1102  /* sw = 3 */
1103  if (sst >= sw3lo && sst < sw2lo) {
1104  weit = (sst - sw3lo) / dtsw3;
1105  pktran[0] = 3;
1106  pktran[1] = 2;
1107  }
1108 
1109  /* sw = 4 */
1110  if (sst < sw3lo) {
1111  weit = 1.;
1112  pktran[0] = 3;
1113  pktran[1] = 3;
1114  }
1115 
1116 
1117  /* Loop through twice, to blend packaging effects */
1118  for (it = 0; it <= 1; it++) {
1119 
1120  pk = pktran[it];
1121 
1122  switch (pk) {
1123  case 1: /* unpackaged */
1124  c0 = upc0;
1125  c1 = upc1;
1126  c2 = upc2;
1127  c3 = upc3;
1128  break;
1129  case 2: /* packaged */
1130  c0 = pc0;
1131  c1 = pc1;
1132  c2 = pc2;
1133  c3 = pc3;
1134  break;
1135  case 3: /*hipackaged */
1136  c0 = hpc0;
1137  c1 = hpc1;
1138  c2 = hpc2;
1139  c3 = hpc3;
1140  break;
1141  default: /* global */
1142  c0 = gc0;
1143  c1 = gc1;
1144  c2 = gc2;
1145  c3 = gc3;
1146  break;
1147  }
1148 
1149 
1150  /* calculate default chl, note coeffs depend on packaging */
1151  exponent = MAX(MIN(c0 + c1 * abr35 + c2 * abr35 * abr35 + c3 * abr35 * abr35*abr35, 10.0), -10.0);
1152  tchl = pow(10.0, exponent);
1153 
1154  /* weight results */
1155  if (it == 0) {
1156  chl = tchl;
1157  } else {
1158  chl = chl * (1. - weit) + tchl * weit;
1159  }
1160  }
1161 
1162  *chlPtr = chl;
1163 
1164  return 0;
1165 
1166 }
1167 
1168 void chl_carder_empirical(l2str *l2rec, float prod[]) {
1169  static int *indx;
1170  static int firstCall = 1;
1171 
1172  int ip, ic;
1173  int32_t ipb;
1174  float ndt;
1175  float sst;
1176  int status;
1177  float rrs[6];
1178  float chl;
1179 
1180  l1str *l1rec = l2rec->l1rec;
1181  int32_t nbands = l1rec->l1file->nbands;
1182 
1183  if (firstCall) {
1184  firstCall = 0;
1185  if ((indx = (int *) calloc(nbands, sizeof (int))) == NULL) {
1186  printf("-E- : Error allocating memory to indx\n");
1187  exit(FATAL_ERROR);
1188  }
1189  indx[0] = bindex_get(412);
1190  indx[1] = bindex_get(443);
1191  indx[2] = bindex_get(490);
1192  indx[3] = bindex_get(510);
1193  if (indx[3] < 0) indx[3] = bindex_get(520);
1194  if (indx[3] < 0) indx[3] = bindex_get(530);
1195  indx[4] = bindex_get(551);
1196  indx[5] = bindex_get(670);
1197  }
1198 
1199  for (ip = 0; ip < l1rec->npix; ip++) {
1200 
1201  ipb = ip*nbands;
1202 
1203  /* grab the remote-sensing reflectances for this pixel */
1204  for (ic = 0; ic < 6; ic++)
1205  rrs[ic] = (l2rec->Rrs[ipb + indx[ic]]);
1206 
1207  /* Use retrieved sst, OR reference sst if no retrieval */
1208  if (l2rec->sst != NULL && l2rec->sst[ip] > -2.0)
1209  sst = l2rec->sst[ip];
1210  else
1211  sst = l1rec->sstref[ip];
1212 
1213  /* grab ndt for this lon/lat */
1214  ndt = get_ndt(l1rec->lon[ip], l1rec->lat[ip]);
1215 
1216  /* run the model */
1217  status = carder_empirical(rrs, sst, ndt, &chl);
1218 
1219  /* store results of this pixel into static global arrays */
1220  if (status == 0) {
1221  prod[ip] = chl;
1222  } else {
1223  prod[ip] = BAD_FLT;
1224  l1rec->flags[ip] |= PRODFAIL;
1225  }
1226  }
1227 }
1228 
1229 
#define MAX(A, B)
Definition: swl0_utils.h:26
integer, parameter int16
Definition: cubeio.f90:3
#define MIN(x, y)
Definition: rice.h:169
#define CFLAG_LO555
Definition: carder.c:56
int j
Definition: decode_rs.h:73
int status
Definition: l1_czcs_hdf.c:32
void iops_carder(l2str *l2rec)
Definition: carder.c:1013
float f1(float x)
int16 * get_flags_carder(l2str *l2rec)
Definition: carder.c:948
#define NX
Definition: carder.c:40
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
#define NDTNX
#define NULL
Definition: decode_rs.h:63
#define CAT_bb_carder
Definition: l2prod.h:85
read l1rec
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
float f3(float z)
#define N_ITER
Definition: carder.c:41
float * lat
#define PRODFAIL
Definition: l2_flags.h:41
#define CFLAG_FAIL
Definition: carder.c:46
#define CAT_aph_carder
Definition: l2prod.h:86
void get_carder(l2str *l2rec, l2prodstr *p, float prod[])
Definition: carder.c:957
int carder_model(int32_t nbands, double *rrs, double *rrswl, float sst, float ndt, double *atot, double *aph, double *adg, double *bb, double *bbp, double *aphmPtr, double *adgmPtr, double *chlPtr, int16 *flags)
Definition: carder.c:100
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed file
Definition: HISTORY.txt:413
#define CAT_chl_carder
Definition: l2prod.h:88
#define CFLAG_DEFAULT
Definition: carder.c:47
#define CFLAG_LO412
Definition: carder.c:55
float f2(float y)
#define CFLAG_BLEND
Definition: carder.c:48
int bindex_get(int32_t wave)
Definition: windex.c:45
#define CAT_bbp_carder
Definition: l2prod.h:99
read recnum
#define CFLAG_PKG
Definition: carder.c:51
#define CAT_adg_carder
Definition: l2prod.h:87
void chl_carder_empirical(l2str *l2rec, float prod[])
Definition: carder.c:1168
float32 slope[]
Definition: l2lists.h:30
float aw_spectra(int32_t wl, int32_t width)
#define FATAL_ERROR
Definition: swl0_parms.h:5
float bbw_spectra(int32_t wl, int32_t width)
const double delta
#define N_PRM_HEADERS
Definition: carder.c:39
integer, parameter double
#define NDTNY
#define CFLAG_GLOBAL
Definition: carder.c:53
int carder_ran(int recnum)
Definition: carder.c:798
#define BAD_FLT
Definition: jplaeriallib.h:19
flags
Definition: DDAlgorithm.h:22
#define CFLAG_UNPKG
Definition: carder.c:50
int32_t nbands
int carder_empirical(float *rrs, float sst, float ndt, float *chlPtr)
Definition: carder.c:1034
#define fabs(a)
Definition: misc.h:93
#define CFLAG_HIPKG
Definition: carder.c:52
Extra metadata that will be written to the HDF4 file l2prod rank
#define CFLAG_CHLINC
Definition: carder.c:57
#define CARDER_FAIL
Definition: carder.c:43
float * lon
data_t s[NROOTS]
Definition: decode_rs.h:75
#define BANDW
Definition: l1.h:52
#define APH_MIN
Definition: carder.c:42
l2prod offset
void carder(l2str *l2rec)
Definition: carder.c:809
int i
Definition: decode_rs.h:71
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int npix
Definition: get_cmp.c:27
#define CAT_a_carder
Definition: l2prod.h:84
float p[MODELMAX]
Definition: atrem_corl1.h:131
void alloc_carder(int32_t npix, int32_t nbands)
Definition: carder.c:784