ocssw V2020
brdf.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 #include "freearray.h"
3 
4 static float radeg = RADEG;
5 static float nw = 1.334;
6 
7 void foq_morel(int foqopt, l2str *l2rec, float wave[], int32_t nwave, float chl,
8  float nLw[], float Fo[], float solz, float senz, float phi, float brdf[]);
9 void dtran_brdf(l2str *l2rec, int32_t ip, float wave[], int32_t nwave, float Fo[], float nLw[], float chl,
10  float brdf[]);
11 void diff_tran_corr_(int *iphase, float *solz, float *senz, float *phi,
12  float *chl, float *taua, float *correct);
13 int32_t getncDimensionLength(int ncid, int dimId);
14 int getncFQdim(int ncid, char *file, char *sdsname, int sds_id, int nexp, float *bdat);
15 
16 /* ---------------------------------------------------------------------------- */
17 /* ocbrdf - bidirectional reflectance correction for one ocean pixel */
18 /* */
19 /* Input: */
20 /* sensorID - sensor identification (for chlorophyll algorithm) */
21 /* brdf_opt - brdf option (bitmask, 1=fresnel(sen), 2=fresnel(sol), 4=f/Q) */
22 /* wave - list of sensor wavelengths (no missing-band gaps) */
23 /* nwave - number of sensor wavelengths to compute brdf */
24 /* solz - solar zenith (deg) */
25 /* senz - sensor zenith (deg) */
26 /* phi - relative azimuth (deg, Gordon definition) */
27 /* nLw - list of normalized water-leaving radiances, band indexed */
28 /* Fo - list of solar irradinaces, band indexed (mW/cm^2/um) */
29 /* */
30 /* Output: */
31 /* brdf - list of brdf corrections to multiply nLw, band indexed */
32 /* */
33 /* Note: for brdf_opt=1, 2, or 3, nLw and Fo can be set to NULL pointer */
34 /* */
35 /* Written By: */
36 /* B. Franz, OBPG, October 2004 */
37 
38 /* ---------------------------------------------------------------------------- */
39 
40 int ocbrdf(l2str *l2rec, int32_t ip, int32_t brdf_opt, float wave[], int32_t nwave,
41  float solz, float senz, float phi, float ws, float chl, float nLw[], float Fo[], float brdf[]) {
42  static int firstCall = 1;
43 
44  float *temp;
45  int32_t status = 0;
46  int32_t iw;
47 
48  if (firstCall == 1) {
49  firstCall = 0;
50  if (brdf_opt > 0)
51  printf("\nApplying ocean BRDF including:\n");
52  if ((brdf_opt & FRESNSEN) > 0)
53  printf(" Reflection/refraction for upwelling radiance.\n");
54  if ((brdf_opt & FRESNSOL) > 0)
55  printf(" Reflection/refraction for downwelling radiance.\n");
56  if ((brdf_opt & FOQMOREL) > 0)
57  printf(" Morel f/Q\n");
58  if ((brdf_opt & QMOREL) > 0)
59  printf(" Morel Q\n");
60  if ((brdf_opt & DTBRDF) > 0)
61  printf(" Gordon diffuse transmittance effect.\n");
62  }
63 
64  if ((temp = (float *) calloc(nwave, sizeof (float))) == NULL) {
65  printf("-E- : Error allocating memory to temp\n");
66  exit(FATAL_ERROR);
67  }
68 
69  for (iw = 0; iw < nwave; iw++)
70  brdf[iw] = 1.0;
71 
72  /* transmittance of view path through air & sea interface */
73  if ((brdf_opt & FRESNSEN) > 0) {
74  float tf = fresnel_sen(senz, 0);
75  for (iw = 0; iw < nwave; iw++)
76  brdf[iw] *= tf;
77  }
78 
79  /* transmittance of solar path through air & sea interface */
80  if ((brdf_opt & FRESNSOL) > 0) {
81  fresnel_sol(wave, nwave, solz, ws, temp, 0);
82  for (iw = 0; iw < nwave; iw++)
83  brdf[iw] *= temp[iw];
84  }
85 
86  /* Morel f/Q correction */
87  if ((brdf_opt & FOQMOREL) > 0) {
88  foq_morel(FOQMOREL, l2rec, wave, nwave, chl, nLw, Fo, solz, senz, phi, temp);
89  for (iw = 0; iw < nwave; iw++) {
90  brdf[iw] *= temp[iw];
91  }
92  }
93 
94  /* Morel f/Q with free sun (Voss) */
95  if ((brdf_opt & QMOREL) > 0) {
96  foq_morel(QMOREL, l2rec, wave, nwave, chl, nLw, Fo, solz, senz, phi, temp);
97  for (iw = 0; iw < nwave; iw++) {
98  brdf[iw] *= temp[iw];
99  }
100  }
101 
102  /* Gordon correction of diffuse transmittance */
103  if ((brdf_opt & DTBRDF) > 0) {
104  dtran_brdf(l2rec, ip, wave, nwave, Fo, nLw, chl, temp);
105  for (iw = 0; iw < nwave; iw++) {
106  brdf[iw] *= temp[iw];
107  }
108  }
109 
110  free(temp);
111 
112  return (status);
113 }
114 
115 /* ---------------------------------------------------------------------------- */
116 /* fresnel_sen() - effects of the air-sea transmittance for sensor view */
117 /* */
118 /* Description: */
119 /* This computes effects of the air-sea transmittance (depending on sensor */
120 /* zenith angle) on the derived normalized water-leaving radiance. */
121 /* Menghua Wang 5/27/02. */
122 /* */
123 /* modified to return fresnel transmittance as option, December 2008, BAF */
124 
125 /* ---------------------------------------------------------------------------- */
126 
127 float fresnel_sen(float senz, int return_tf) {
128  static float tf0 = 0.9795218;
129 
130  float mu = cos(senz / radeg);
131  float sq, r2, q1, fres, tf, brdf;
132 
133  sq = sqrt(nw * nw - 1. + mu * mu);
134  r2 = pow((mu - sq) / (mu + sq), 2);
135  q1 = (1. - mu * mu - mu * sq) / (1. - mu * mu + mu * sq);
136  fres = r2 * (q1 * q1 + 1.) / 2.0;
137  tf = 1. - fres;
138  brdf = tf0 / tf;
139 
140  if (return_tf != 0)
141  return (tf);
142  else
143  return (brdf);
144 }
145 
146 
147 /* ---------------------------------------------------------------------------- */
148 /* fresnel_sol() - effects of the air-sea transmittance for solar path */
149 /* */
150 /* Description: */
151 /* This computes the correction factor on normalized water-leaving radiance */
152 /* to account for the solar zenith angle effects on the downward irradiance */
153 /* from above the ocean surface to below the surface. */
154 /* Menghua Wang 9/27/04. */
155 /* */
156 /* Added windspeed dependence, December 2004, BAF */
157 /* Modified to return air-sea transmittance as option, December 2008, BAF */
158 
159 /* ---------------------------------------------------------------------------- */
160 void fresnel_sol(float wave[], int32_t nwave, float solz, float ws, float brdf[], int return_tf) {
161 
162 #define NTWAVE 6
163 #define NTWIND 5
164 
165  static int firstCall = 1;
166  static int *tindx;
167  static float *tf0;
168  static float twave [] = {412., 443., 490., 510., 555., 670.};
169  static float tsigma[] = {0.0, 0.1, 0.2, 0.3, 0.4};
170 
171  /* M Wang, personal communication, red-nir iterpolated */
172  static float tf0_w[] = {412., 443., 490., 510., 555., 670., 765., 865.};
173  static float tf0_v[] = {0.965980, 0.968320, 0.971040, 0.971860, 0.973450, 0.977513, 0.980870, 0.984403};
174 
175  static float c[NTWIND][4][NTWAVE] = {
176  { /* ws=0.0 */
177  { -0.0087, -0.0122, -0.0156, -0.0163, -0.0172, -0.0172}, /* a */
178  { 0.0638, 0.0415, 0.0188, 0.0133, 0.0048, -0.0003}, /* b */
179  { -0.0379, -0.0780, -0.1156, -0.1244, -0.1368, -0.1430}, /* c */
180  { -0.0311, -0.0427, -0.0511, -0.0523, -0.0526, -0.0478} /* d */
181  },
182  { /* ws=1.9 */
183  { -0.0011, -0.0037, -0.0068, -0.0077, -0.0090, -0.0106}, /* a */
184  { 0.0926, 0.0746, 0.0534, 0.0473, 0.0368, 0.0237}, /* b */
185  { -5.3E-4, -0.0371, -0.0762, -0.0869, -0.1048, -0.1260}, /* c */
186  { -0.0205, -0.0325, -0.0438, -0.0465, -0.0506, -0.0541} /* d */
187  },
188  { /* ws=7.5 */
189  { 6.8E-5, -0.0018, -0.0011, -0.0012, -0.0015, -0.0013}, /* a */
190  { 0.1150, 0.1115, 0.1075, 0.1064, 0.1044, 0.1029}, /* b */
191  { 0.0649, 0.0379, 0.0342, 0.0301, 0.0232, 0.0158}, /* c */
192  { 0.0065, -0.0039, -0.0036, -0.0047, -0.0062, -0.0072} /* d */
193  },
194  { /* ws=16.9 */
195  { -0.0088, -0.0097, -0.0104, -0.0106, -0.0110, -0.0111}, /* a */
196  { 0.0697, 0.0678, 0.0657, 0.0651, 0.0640, 0.0637}, /* b */
197  { 0.0424, 0.0328, 0.0233, 0.0208, 0.0166, 0.0125}, /* c */
198  { 0.0047, 0.0013, -0.0016, -0.0022, -0.0031, -0.0036} /* d */
199  },
200  { /* ws=30 */
201  { -0.0081, -0.0089, -0.0096, -0.0098, -0.0101, -0.0104}, /* a */
202  { 0.0482, 0.0466, 0.0450, 0.0444, 0.0439, 0.0434}, /* b */
203  { 0.0290, 0.0220, 0.0150, 0.0131, 0.0103, 0.0070}, /* c */
204  { 0.0029, 0.0004, -0.0017, -0.0022, -0.0029, -0.0033} /* d */
205  }
206  };
207 
208  float x, x2, x3, x4;
209  int is, is1, is2;
210  int iw, i;
211  float sigma;
212  float slp;
213  float brdf1, brdf2;
214 
215  /* on first call, find closest table entry to each input wavelength */
216  if (firstCall == 1) {
217  firstCall = 0;
218  if ((tindx = (int *) calloc(nwave, sizeof (int))) == NULL) {
219  printf("-E- : Error allocating memory to tindx\n");
220  exit(FATAL_ERROR);
221  }
222  if ((tf0 = (float *) calloc(nwave, sizeof (float))) == NULL) {
223  printf("-E- : Error allocating memory to tf0\n");
224  exit(FATAL_ERROR);
225  }
226 
227  for (iw = 0; iw < nwave; iw++) {
228  tindx[iw] = windex(wave[iw], twave, NTWAVE);
229  tf0 [iw] = linterp(tf0_w, tf0_v, 8, wave[iw]);
230  }
231  }
232 
233  sigma = 0.0731 * sqrt(MAX(ws, 0.0));
234 
235  x = log(cos(MIN(solz, 80.0) / radeg));
236  x2 = x*x;
237  x3 = x*x2;
238  x4 = x*x3;
239 
240  /* find bracketing table winds */
241  for (is = 0; is < NTWIND; is++)
242  if (tsigma[is] > sigma)
243  break;
244  is2 = MIN(is, NTWIND - 1);
245  is1 = is2 - 1;
246  slp = (sigma - tsigma[is1]) / (tsigma[is2] - tsigma[is1]);
247 
248  /* compute at bounding windspeeds and interpolate */
249  for (iw = 0; iw < nwave; iw++) {
250  i = tindx[iw];
251  brdf1 = 1.0 + c[is1][0][i] * x + c[is1][1][i] * x2
252  + c[is1][2][i] * x3 + c[is1][3][i] * x4;
253  brdf2 = 1.0 + c[is2][0][i] * x + c[is2][1][i] * x2
254  + c[is2][2][i] * x3 + c[is2][3][i] * x4;
255  brdf[iw] = brdf1 + slp * (brdf2 - brdf1);
256  if (return_tf != 0)
257  brdf[iw] = tf0[iw] / brdf[iw];
258  }
259 }
260 
261 
262 /* ---------------------------------------------------------------------------- */
263 /* gothic_R - returns total effect of transmission through air-sea interface */
264 
265 /* ---------------------------------------------------------------------------- */
266 void gothic_R(float wave[], int32_t nwave, float solz, float senz, float ws, float R[]) {
267  int iw;
268  float tf = fresnel_sen(senz, 1);
269 
270  fresnel_sol(wave, nwave, solz, ws, R, 1);
271  for (iw = 0; iw < nwave; iw++)
272  R[iw] = R[iw] * tf / nw / nw;
273 
274  return;
275 }
276 
277 
278 #define N_W 7
279 #define N_S 6
280 #define N_C 6
281 #define N_N 17
282 #define N_A 13
283 
284 /* return closest indice of xtab where xtab[i] < x */
285 int morel_index(float xtab[], int32_t ntab, float x) {
286  int i = 0;
287 
288  if (x <= xtab[0])
289  i = 0;
290  else if (x >= xtab[ntab - 1])
291  i = ntab - 2;
292  else {
293  while (x > xtab[i]) i++;
294  i--;
295  }
296 
297  return (i);
298 }
299 
300 
301 /* ---------------------------------------------------------------------------- */
302 /* foqint_morel() - reads and interpolates f/Q tables of Morel & Gentilli */
303 /* */
304 /* wave[] - list of input wavelengths (nm) */
305 /* nwave - number of input wavelengths */
306 /* solz - solar zenith angle of observation (deg) */
307 /* senzp - view zenith angle of observation, below surface (deg) */
308 /* phi - relative azimuth of observation, 0=<--> (deg) */
309 /* chl - chlorophyll-a concentration (mg/m^3) */
310 /* brdf[] - band-indexed array of f/Q corrections per wavelength (f0/Q0)/(f/Q) */
311 /* */
312 
313 /* ---------------------------------------------------------------------------- */
314 void foqint_morel(char *file, float wave[], int32_t nwave, float solz, float senzp,
315  float phi, float chl, float brdf[]) {
316  static int firstCall = 1;
317  static float *foqtab;
318  static float *wavetab;
319  static float *solztab;
320  static float *chltab;
321  static float *senztab;
322  static float *phitab;
323  static float *lchltab;
324 
325 
326  float lchl;
327  int i, j, k, l, m;
328  int iw, js, kc, ln, ma;
329  float ds[2], dc[2], dn[2], da[2];
330 
331  static int32_t n_a, n_n, n_c, n_s, n_w;
332 
333  char sdsname[H4_MAX_NC_NAME];
334  int ncid, ndims, nvars, ngatts, unlimdimid;
335  int32 sds_id;
336  nc_type rh_type; /* variable type */
337  int rh_dimids[H4_MAX_VAR_DIMS]; /* dimension IDs */
338  int rh_natts; /* number of attributes */
339  int rh_ndims; /* number of dims */
340  int status, ndx;
341  if (firstCall == 1) {
342 
343  /* try netCDF first */
344  if (nc_open(file, NC_NOWRITE, &ncid) == NC_NOERR) {
345 
346  status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid);
347  if (status != NC_NOERR) {
348  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
349  __FILE__, __LINE__, sdsname, file);
350  exit(1);
351  }
352 
353  strcpy(sdsname, "foq");
354 
355  status = nc_inq_varid(ncid, sdsname, &sds_id);
356  if (status != NC_NOERR) {
357  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
358  __FILE__, __LINE__, sdsname, file);
359  exit(1);
360  }
361  status = nc_inq_var(ncid, sds_id, 0, &rh_type, &rh_ndims, rh_dimids,
362  &rh_natts);
363 
364  if (rh_ndims != 5) {
365  fprintf(stderr, "-E- %s line %d: Wrong number of dimensions for %s. Need 5 got %d.\n",
366  __FILE__, __LINE__, sdsname, rh_ndims);
367  exit(1);
368 
369  }
370  printf("rh_ndims=%d rh_dimids=%d %d %d %d %d \n", rh_ndims, rh_dimids[0], rh_dimids[1], rh_dimids[2], rh_dimids[3], rh_dimids[4]);
371 
372  n_w = getncDimensionLength(ncid, rh_dimids[0]);
373  n_s = getncDimensionLength(ncid, rh_dimids[1]);
374  n_c = getncDimensionLength(ncid, rh_dimids[2]);
375  n_n = getncDimensionLength(ncid, rh_dimids[3]);
376  n_a = getncDimensionLength(ncid, rh_dimids[4]);
377 
378  printf("morel f/q file dimensions n_a=%d n_n=%d n_c=%d n_s=%d n_w=%d \n", n_a, n_n, n_c, n_s, n_w);
379 
380  printf("\nReading foq file %s ndims=%d nvars=%d sds_id=%d var=%s\n\n", file, ndims, nvars, sds_id, sdsname);
381  /* Read the data. */
382  if ((foqtab = (float *) calloc(n_w * n_s * n_c * n_n * n_a, sizeof (float))) == NULL) {
383  printf("-E- : Error allocating memory to tindx\n");
384  exit(FATAL_ERROR);
385  }
386  if (nc_get_var(ncid, sds_id, foqtab) != NC_NOERR) {
387  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
388  __FILE__, __LINE__, sdsname, file);
389  exit(1);
390  }
391 
392  if ((phitab = (float *) calloc(n_a, sizeof (float))) == NULL) {
393  printf("-E- : Error allocating memory to tindx\n");
394  exit(FATAL_ERROR);
395  }
396  if ((senztab = (float *) calloc(n_n, sizeof (float))) == NULL) {
397  printf("-E- : Error allocating memory to tindx\n");
398  exit(FATAL_ERROR);
399  }
400  if ((solztab = (float *) calloc(n_s, sizeof (float))) == NULL) {
401  printf("-E- : Error allocating memory to tindx\n");
402  exit(FATAL_ERROR);
403  }
404  if ((chltab = (float *) calloc(n_c, sizeof (float))) == NULL) {
405  printf("-E- : Error allocating memory to tindx\n");
406  exit(FATAL_ERROR);
407  }
408 
409  if ((wavetab = (float *) calloc(n_w, sizeof (float))) == NULL) {
410  printf("-E- : Error allocating memory to tindx\n");
411  exit(FATAL_ERROR);
412  }
413  if ((lchltab = (float *) calloc(n_c, sizeof (float))) == NULL) {
414  printf("-E- : Error allocating memory to tindx\n");
415  exit(FATAL_ERROR);
416  }
417 
418  strcpy(sdsname, "phi");
419  status = getncFQdim(ncid, file, sdsname, sds_id, n_a, phitab);
420  strcpy(sdsname, "senz");
421  status = getncFQdim(ncid, file, sdsname, sds_id, n_n, senztab);
422  strcpy(sdsname, "solz");
423  status = getncFQdim(ncid, file, sdsname, sds_id, n_s, solztab);
424  strcpy(sdsname, "chl");
425  status = getncFQdim(ncid, file, sdsname, sds_id, n_c, chltab);
426  strcpy(sdsname, "wave");
427  status = getncFQdim(ncid, file, sdsname, sds_id, n_w, wavetab);
428 
429  printf("\nClosing foq file %s\n\n", file);
430 
431  /* Close the file */
432  if (nc_close(ncid) != NC_NOERR) {
433  fprintf(stderr, "-E- %s line %d: error closing %s.\n",
434  __FILE__, __LINE__, file);
435  exit(1);
436  }
437  }
438 
439 
440  printf("\nMorel f/Q table from file %s\n", file);
441 
442  // for (i=0; i<n_w; i++)
443  // for (j=0; j<n_s; j++)
444  // for (k=0; k<n_c; k++)
445  // for (l=0; l<n_n; l++) {
446  // for (m=0; m<n_a-1; m++)
447  // {
448  // ndx = i*n_s*n_c*n_n*n_a+j*n_c*n_n*n_a+k*n_n*n_a+l*n_a+m;
449  // printf("%f",*(foqtab+ndx));
450  // //printf("%f",foqtab[i][j][k][l][m]);
451  // }
452  // m = n_a-1;
453  // ndx = i*n_s*n_c*n_n*n_a+j*n_c*n_n*n_a+k*n_n*n_a+l*n_a+m;
454  // //printf("%f\n",foqtab[i][j][k][l][m]);
455  // printf("%f\n",*(foqtab+ndx));
456  // }
457 
458  /* create table for log of chlorophyll */
459  for (kc = 0; kc < n_c; kc++)
460  lchltab[kc] = log(chltab[kc]);
461 
462  firstCall = 0;
463  }
464 
465 
466  lchl = log(MAX(chl, 0.01));
467 
468  if (senzp < senztab[0]) {
469  senzp = senztab[0];
470  }
471 
472  /* lower bounding indices */
473  js = morel_index(solztab, n_s, solz);
474  kc = morel_index(lchltab, n_c, lchl);
475  ln = morel_index(senztab, n_n, senzp);
476  ma = morel_index(phitab, n_a, phi);
477 
478  ds[0] = (solztab[js + 1] - solz) / (solztab[js + 1] - solztab[js]);
479  ds[1] = (solz - solztab[js]) / (solztab[js + 1] - solztab[js]);
480 
481  dc[0] = (lchltab[kc + 1] - lchl) / (lchltab[kc + 1] - lchltab[kc]);
482  dc[1] = (lchl - lchltab[kc]) / (lchltab[kc + 1] - lchltab[kc]);
483 
484  dn[0] = (senztab[ln + 1] - senzp) / (senztab[ln + 1] - senztab[ln]);
485  dn[1] = (senzp - senztab[ln]) / (senztab[ln + 1] - senztab[ln]);
486 
487  da[0] = (phitab [ma + 1] - phi) / (phitab [ma + 1] - phitab [ma]);
488  da[1] = (phi - phitab [ma]) / (phitab [ma + 1] - phitab [ma]);
489 
490  for (iw = 0; iw < nwave; iw++) {
491 
492  /* using nearest wavelength (tables are for MERIS bands) */
493 
494  i = windex(wave[iw], wavetab, n_w);
495 
496  brdf[iw] = 0.0;
497 
498  for (j = 0; j <= 1; j++)
499  for (k = 0; k <= 1; k++)
500  for (l = 0; l <= 1; l++)
501  for (m = 0; m <= 1; m++) {
502 
503  ndx = i * n_s * n_c * n_n * n_a + (js + j) * n_c * n_n * n_a + (kc + k) * n_n * n_a + (ln + l) * n_a + ma + m;
504  brdf[iw] += ds[j] * dc[k] * dn[l] * da[m]*(*(foqtab + ndx));
505  }
506  }
507 
508  return;
509 }
510 
511 /* -----------------------------------------------------------------------*/
512 /* getDimensionLength() - get the length of the dimension ID */
513 
514 /* -----------------------------------------------------------------------*/
515 int32_t getncDimensionLength(int ncid, int dimId) {
516  size_t length;
517  int32_t nl;
518 
519  if (nc_inq_dimlen(ncid, dimId, &length) != NC_NOERR) {
520  char name[H4_MAX_NC_NAME];
521  nc_inq_dim(ncid, dimId, name, &length);
522  fprintf(stderr,
523  "-E- %s line %d: could not get size of dimension \"%s\" in netCDF File.\n",
524  __FILE__, __LINE__, name);
525  exit(1);
526  }
527  nl = length;
528  return nl;
529 }
530 
531 int getncFQdim(int ncid, char *file, char *sdsname, int sds_id, int nexp, float *bdat) {
532 
533  int status;
534  nc_type rh_type; /* variable type */
535  int rh_dimids[H4_MAX_VAR_DIMS]; /* dimension IDs */
536  int rh_natts; /* number of attributes */
537  int rh_ndims, nsz; /* number of dims */
538 
539  status = nc_inq_varid(ncid, sdsname, &sds_id);
540  if (status != NC_NOERR) {
541  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
542  __FILE__, __LINE__, sdsname, file);
543  exit(1);
544  }
545  status = nc_inq_var(ncid, sds_id, 0, &rh_type, &rh_ndims, rh_dimids,
546  &rh_natts);
547 
548  nsz = getncDimensionLength(ncid, rh_dimids[0]);
549 
550  if (nsz != nexp) {
551  fprintf(stderr, "-E- %s line %d: Wrong dimensions for %s. Need %d got %d.\n",
552  __FILE__, __LINE__, sdsname, nexp, nsz);
553  exit(1);
554 
555 
556  }
557  if (nc_get_var(ncid, sds_id, bdat) != NC_NOERR) {
558  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
559  __FILE__, __LINE__, sdsname, file);
560  exit(1);
561  }
562 
563  return status;
564 
565 }
566 /* ---------------------------------------------------------------------------- */
567 /* foq_morel() - computes f/Q correction of Morel & Gentilli by iteration */
568 /* */
569 /* foqopt - 0=full f/Q, 1=no normalization to sun overhead (fixed f) */
570 /* sensorID - MSl12 sensor identification number */
571 /* wave[] - list of input wavelengths (nm) */
572 /* nwave - number of input wavelengths */
573 /* nLw[] - normalize water-leaving radiances per wave (mW/cm^2/um/sr) */
574 /* Fo[] - solar irradiance per wave (mW/cm^2/um/sr) */
575 /* solz - solar zenith angle of observation (deg) */
576 /* senzp - view zenith angle of observation, below surface (deg) */
577 /* phi - relative azimuth of observation, 0=<--> (deg) */
578 /* brdf[] - band-indexed array of f/Q corrections per wavelength (f0/Q0)/(f/Q) */
579 
580 /* ---------------------------------------------------------------------------- */
581 void foq_morel(int foqopt, l2str *l2rec, float wave[], int32_t nwave, float chl,
582  float nLw[], float Fo[], float solz, float senz, float phi, float brdf[]) {
583  static int maxiter = 3;
584  static int compchl = 1;
585 
586  float senzp, phip;
587  float *foq0;
588  float *foq;
589  float *Rrs;
590  int32_t iw, iter;
591  int numiter;
592 
593  if ((foq0 = (float *) calloc(nwave, sizeof (float))) == NULL) {
594  printf("-E- : Error allocating memory to foq0\n");
595  exit(FATAL_ERROR);
596  }
597  if ((foq = (float *) calloc(nwave, sizeof (float))) == NULL) {
598  printf("-E- : Error allocating memory to foq\n");
599  exit(FATAL_ERROR);
600  }
601  if ((Rrs = (float *) calloc(nwave, sizeof (float))) == NULL) {
602  printf("-E- : Error allocating memory to Rrs\n");
603  exit(FATAL_ERROR);
604  }
605 
606  /* Need view zenith and relative azimuth below water. The MSl12 definition of */
607  /* relative azimuth is consistent with the below-water definition. We just need */
608  /* to limit from [-180,180] to [0,180]. */
609  phip = fabs(phi);
610  senzp = asin(sin(senz / radeg) / nw) * radeg;
611 
612  /* Compute starting chlorophyll (if not supplied) */
613 
614  if (chl < 0.0) {
615  compchl = 1;
616  numiter = maxiter;
617  for (iw = 0; iw < nwave; iw++) {
618  Rrs[iw] = nLw[iw] / Fo[iw];
619  }
620  chl = get_default_chl(l2rec, Rrs);
621  } else {
622  compchl = 0;
623  numiter = 1;
624  }
625 
626  /* If we retrieved a valid chlorophyll, we can compute the correction. */
627 
628  if (chl >= 0.0) {
629 
630  for (iter = 0; iter < numiter; iter++) {
631 
632  if (foqopt == QMOREL)
633  foqint_morel(input->fqfile, wave, nwave, solz, 0.0, 0.0, chl, foq0);
634  else
635  foqint_morel(input->fqfile, wave, nwave, 0.0, 0.0, 0.0, chl, foq0);
636 
637  foqint_morel(input->fqfile, wave, nwave, solz, senzp, phip, chl, foq);
638 
639  for (iw = 0; iw < nwave; iw++) {
640  brdf [iw] = foq0[iw] / foq[iw];
641  Rrs [iw] = nLw[iw] * brdf[iw] / Fo[iw];
642  }
643 
644  if (compchl) {
645  chl = get_default_chl(l2rec, Rrs);
646  }
647  }
648  }
649 
650  if (chl < 0.0) {
651  for (iw = 0; iw < nwave; iw++) {
652  brdf [iw] = 1.0;
653  }
654  }
655 
656  free(foq);
657  free(foq0);
658  free(Rrs);
659  return;
660 }
661 
662 
663 /* ---------------------------------------------------------------------------- */
664 /* qint_morel() - reads and interpolates Q tables of Morel & Gentilli */
665 /* */
666 /* wave[] - list of input wavelengths (nm) */
667 /* nwave - number of input wavelengths */
668 /* solz - solar zenith angle of observation (deg) */
669 /* */
670 
671 /* ---------------------------------------------------------------------------- */
672 void qint_morel(float wave[], int32_t nwave, float solz, float chl, float Qn[]) {
673  static int firstCall = 1;
674 
675  static float q0tab [N_C][N_W];
676  static float sqtab [N_C][N_W];
677  static float dq0tab [N_C][N_W];
678  static float dsqtab [N_C][N_W];
679  static float wavetab[N_W] = {412.5, 442.5, 490.0, 510.0, 560.0, 620.0, 660.0};
680  static float chltab [N_C] = {0.03, 0.1, 0.3, 1.0, 3.0, 10.0};
681  static float lchltab[N_C];
682 
683  int ic, iw;
684  float lchl, q1, q2, s1, s2, qn1, qn2;
685 
686  if (firstCall == 1) {
687 
688  char filename[FILENAME_MAX];
689  char *delim = " ";
690  FILE *fp;
691  char *tmp_str;
692  char line [80];
693  char buff [80];
694 
695  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
696  printf("OCDATAROOT environment variable is not defined.\n");
697  exit(1);
698  }
699 
700  strcpy(filename, tmp_str);
701  strcat(filename, "/common/morel_q0.dat");
702  fp = fopen(filename, "r");
703  if (!fp) {
704  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, filename);
705  exit(1);
706  }
707 
708  printf("\nLoading Morel Q0 table from file %s\n", filename);
709 
710  for (ic = 0; ic < N_C; ic++) {
711  if (fgets(line, 80, fp) == NULL) {
712  fprintf(stderr, "-E- %s line %d: error reading %s at line %d\n", __FILE__, __LINE__, filename, ic);
713  exit(1);
714  }
715  strcpy(buff, line);
716  chltab[ic] = atof(strtok(buff, delim));
717  for (iw = 0; iw < N_W; iw++) {
718  q0tab[ic][iw] = atof(strtok(NULL, delim));
719  }
720  }
721 
722  strcpy(filename, tmp_str);
723  strcat(filename, "/common/morel_sq.dat");
724  fp = fopen(filename, "r");
725  if (!fp) {
726  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, filename);
727  exit(1);
728  }
729 
730  printf("\nLoading Morel Sq table from file %s\n", filename);
731 
732  for (ic = 0; ic < N_C; ic++) {
733  if (fgets(line, 80, fp) == NULL) {
734  fprintf(stderr, "-E- %s line %d: error reading %s at line %d\n", __FILE__, __LINE__, filename, ic);
735  exit(1);
736  }
737  strcpy(buff, line);
738  chltab[ic] = atof(strtok(buff, delim));
739  for (iw = 0; iw < N_W; iw++) {
740  sqtab[ic][iw] = atof(strtok(NULL, delim));
741  }
742  }
743 
744  // create table for log of chlorophyll
745  for (ic = 0; ic < N_C; ic++)
746  lchltab[ic] = log(chltab[ic]);
747 
748  // precompute derivatives for spline interpolation in wavelength
749 
750  for (ic = 0; ic < N_C; ic++) {
751  spline(wavetab, &q0tab[ic][0], N_W, 1e30, 1e30, &dq0tab[ic][0]);
752  spline(wavetab, &sqtab[ic][0], N_W, 1e30, 1e30, &dsqtab[ic][0]);
753  }
754 
755  firstCall = 0;
756  }
757 
758  // find bounding chl indices and weights
759 
760  lchl = log(MAX(chl, 0.01));
761  ic = morel_index(lchltab, N_C, lchl);
762 
763  for (iw = 0; iw < nwave; iw++) {
764 
765  if (wave[iw] <= wavetab[0]) {
766  q1 = q0tab[ic + 0][0];
767  q2 = q0tab[ic + 1][0];
768  s1 = sqtab[ic + 0][0];
769  s2 = sqtab[ic + 1][0];
770  } else if (wave[iw] >= wavetab[N_W - 1]) {
771  q1 = q0tab[ic + 0][N_W - 1];
772  q2 = q0tab[ic + 1][N_W - 1];
773  s1 = sqtab[ic + 0][N_W - 1];
774  s2 = sqtab[ic + 1][N_W - 1];
775  } else {
776  splint(wavetab, &q0tab[ic + 0][0], &dq0tab[ic + 0][0], N_W, wave[iw], &q1);
777  splint(wavetab, &q0tab[ic + 1][0], &dq0tab[ic + 1][0], N_W, wave[iw], &q2);
778  splint(wavetab, &sqtab[ic + 0][0], &dsqtab[ic + 0][0], N_W, wave[iw], &s1);
779  splint(wavetab, &sqtab[ic + 1][0], &dsqtab[ic + 1][0], N_W, wave[iw], &s2);
780  }
781 
782  qn1 = q1 + s1 * (1.0 - cos(solz / radeg));
783  qn2 = q2 + s2 * (1.0 - cos(solz / radeg));
784 
785  Qn[iw] = qn1 + (lchl - lchltab[ic])*(qn2 - qn1) / (lchltab[ic + 1] - lchltab[ic]);
786  }
787 
788  return;
789 }
790 
791 
792 /* ---------------------------------------------------------------------------- */
793 /* dtran_brdf() - computes BRDF effect on diffuse transmittance */
794 /* */
795 /* wrapper for H. Gordon's diff_tran_corr() function. */
796 /* */
797 /* l2rec - pointer to level-2 record (assumes aerosol terms loaded for ip) */
798 /* ip - scan pixel number (from 0) */
799 /* wave[] - list of input wavelengths (nm) */
800 /* nwave - number of input wavelengths */
801 /* Fo[] - solar irradiance per wave (mW/cm^2/um/sr) */
802 /* nLw[] - normalize water-leaving radiances per wave (mW/cm^2/um/sr) */
803 /* brdf[] - band-indexed array of correction to nLw */
804 /* */
805 /* B. Franz, September 2006 */
806 
807 /* ---------------------------------------------------------------------------- */
808 void dtran_brdf(l2str *l2rec, int32_t ip, float wave[], int32_t nwave, float Fo[], float nLw[], float chl,
809  float brdf[]) {
810  static int firstCall = 1;
811  static int modindex[] = {4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15, 16};
812  static float wavetab [] = {412., 443., 490., 510., 555., 670., 765., 865.};
813  static int *waveindex;
814  static int iw865;
815 
816  int32_t amodmin = l2rec->aermodmin[ip];
817  int32_t amodmax = l2rec->aermodmax[ip];
818  float amodrat = l2rec->aerratio[ip];
819  float solz = l2rec->l1rec->solz[ip];
820  float senz = l2rec->l1rec->senz[ip];
821  float phi = l2rec->l1rec->delphi[ip];
822  float taua865;
823 
824  float *Rrs;
825  float **cbrdf;
826  int iphase;
827  float *correct;
828  int iw, im;
829 
830  if (firstCall == 1) {
831  if ((waveindex = (int *) calloc(nwave, sizeof (int))) == NULL) {
832  printf("-E- : Error allocating memory to waveindex\n");
833  exit(FATAL_ERROR);
834  }
835  printf("Initializing BRDF correction to diffuse transmittance.\n");
836  /* tables are for seawifs wavelength; use nearest */
837  for (iw = 0; iw < nwave; iw++) {
838  waveindex[iw] = windex(wave[iw], wavetab, 8);
839  printf("Using %3.0f entries for sensor wavelength %3.0f nm\n",
840  wavetab[waveindex[iw]], wave[iw]);
841  }
842  iw865 = waveindex[7];
843  firstCall = 0;
844  }
845  if ((Rrs = (float *) calloc(nwave, sizeof (float))) == NULL) {
846  printf("-E- : Error allocating memory to Rrs\n");
847  exit(FATAL_ERROR);
848  }
849  if ((correct = (float *) calloc(nwave, sizeof (float))) == NULL) {
850  printf("-E- : Error allocating memory to correct\n");
851  exit(FATAL_ERROR);
852  }
853  if ((cbrdf = (float **) calloc(2, sizeof (float *))) == NULL) {
854  printf("-E- : Error allocating memory to cbrdf\n");
855  exit(FATAL_ERROR);
856  }
857  for (im = 0; im < 2; im++)
858  if ((cbrdf[im] = (float *) calloc(2, sizeof (float))) == NULL) {
859  printf("-E- : Error allocating memory to cbrdf\n");
860  exit(FATAL_ERROR);
861  }
862 
863  /* want relative azimuth in radians, other angles in degrees */
864  phi /= radeg;
865 
866  /* Compute starting chlorophyll (if not supplied) */
867  if (chl < 0.0) {
868  for (iw = 0; iw < nwave; iw++) {
869  Rrs[iw] = nLw[iw] / Fo[iw];
870  }
871  chl = get_default_chl(l2rec, Rrs);
872  }
873 
874  taua865 = l2rec->taua[ip * nwave + iw865];
875 
876  /* If we retrieved a valid chlorophyll, we can compute the correction. */
877 
878  if (chl > 0.0 && taua865 > 0.0) {
879 
880  chl = MAX(MIN(chl, 10.0), 0.03);
881 
882  for (im = 0; im <= 1; im++) {
883 
884  /* 12-model index within 16-model suite (assumes std 12 models) */
885  if (im == 0)
886  iphase = modindex[amodmin];
887  else
888  iphase = modindex[amodmax];
889 
890  diff_tran_corr_(&iphase, &solz, &senz, &phi, &chl, &taua865, correct);
891 
892  /* compute correction per wavelength, per aerosol model */
893  for (iw = 0; iw < nwave; iw++) {
894  cbrdf[im][iw] = (1 + correct[waveindex[iw]]);
895  }
896  }
897 
898  /* combine corrections from 2 aerosol models */
899  for (iw = 0; iw < nwave; iw++) {
900  brdf[iw] = 1.0 / ((1 - amodrat) * cbrdf[0][iw] + amodrat * cbrdf[1][iw]);
901  }
902 
903  } else {
904  for (iw = 0; iw < nwave; iw++)
905  brdf [iw] = 1.0;
906  }
907 
908  free(Rrs);
909  free(correct);
910  freeArrayf(cbrdf, 2);
911  return;
912 }
void gothic_R(float wave[], int32_t nwave, float solz, float senz, float ws, float R[])
Definition: brdf.c:266
#define NTWAVE
#define MAX(A, B)
Definition: swl0_utils.h:26
#define MIN(x, y)
Definition: rice.h:169
void freeArrayf(float **a, int32_t m)
Definition: freearray.h:10
int j
Definition: decode_rs.h:73
void diff_tran_corr_(int *iphase, float *solz, float *senz, float *phi, float *chl, float *taua, float *correct)
int status
Definition: l1_czcs_hdf.c:31
void foq_morel(int foqopt, l2str *l2rec, float wave[], int32_t nwave, float chl, float nLw[], float Fo[], float solz, float senz, float phi, float brdf[])
Definition: brdf.c:581
#define NTWIND
float m
Definition: mie_kernal.py:30
int32_t nl
Definition: atrem_corl1.h:132
constexpr double RADEG
#define FOQMOREL
Definition: l12_parms.h:89
#define NULL
Definition: decode_rs.h:63
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
float mu
int getncFQdim(int ncid, char *file, char *sdsname, int sds_id, int nexp, float *bdat)
Definition: brdf.c:531
int morel_index(float xtab[], int32_t ntab, float x)
Definition: brdf.c:285
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start thereby resolving MODur00108 Changed header guard macro names to avoid reserved name resolving MODur00104 Maneuver list file for Terra satellite was updated to include data from Jecue DuChateu Maneuver list files for both Terra and Aqua were updated to include two maneuvers from recent IOT weekly reports The limits for Z component of angular momentum was and to set the fourth GEO scan quality flag for a scan depending upon whether or not it occurred during one of them Added _FillValue attributes to many and changed the fill value for the sector start times from resolving MODur00072 Writes boundingcoordinate metadata to L1A archived metadata For PERCENT *ECS change to treat rather resolving GSFcd01518 Added a LogReport Changed to create the Average Temperatures vdata even if the number of scans is
Definition: HISTORY.txt:301
#define QMOREL
Definition: l12_parms.h:91
void dtran_brdf(l2str *l2rec, int32_t ip, float wave[], int32_t nwave, float Fo[], float nLw[], float chl, float brdf[])
Definition: brdf.c:808
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
int32_t im
Definition: atrem_corl1.h:161
instr * input
float fresnel_sen(float senz, int return_tf)
Definition: brdf.c:127
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
#define FATAL_ERROR
Definition: swl0_parms.h:5
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
#define N_C
Definition: brdf.c:280
#define N_W
Definition: brdf.c:278
void foqint_morel(char *file, float wave[], int32_t nwave, float solz, float senzp, float phi, float chl, float brdf[])
Definition: brdf.c:314
#define FRESNSOL
Definition: l12_parms.h:88
#define FRESNSEN
Definition: l12_parms.h:87
float get_default_chl(l2str *l2rec, float Rrs[])
Definition: get_chl.c:478
#define fabs(a)
Definition: misc.h:93
char * name
Definition: Granule.c:1234
#define DTBRDF
Definition: l12_parms.h:90
subroutine splint(xa, ya, y2a, n, x, y)
void qint_morel(float wave[], int32_t nwave, float solz, float chl, float Qn[])
Definition: brdf.c:672
void fresnel_sol(float wave[], int32_t nwave, float solz, float ws, float brdf[], int return_tf)
Definition: brdf.c:160
#define R
Definition: make_L3_v1.1.c:96
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 as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT and to use it determine the correct files to retrieve the ephemeris and attitude data from Changed to produce a LocalGranuleID starting with MYD03 if run on Aqua data Added the Scan Type file attribute to the Geolocation copied from the L1A and attitude_angels to radians rather than degrees The accumulation of Cumulated gflags was moved from GEO_validate_earth_location c to GEO_locate_one_scan c
Definition: HISTORY.txt:464
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 ocbrdf(l2str *l2rec, int32_t ip, int32_t brdf_opt, float wave[], int32_t nwave, float solz, float senz, float phi, float ws, float chl, float nLw[], float Fo[], float brdf[])
Definition: brdf.c:40
int32_t getncDimensionLength(int ncid, int dimId)
Definition: brdf.c:515
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:71
int k
Definition: decode_rs.h:73