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
atmocor2.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 #include "atrem_corl1.h"
3 
4 /* --------------------------------------------------------------------------------------- */
5 /* atmocor2() - atmospheric correction, converts Lt -> nLw */
6 /* */
7 /* C-version, September 2004, B. Franz */
8 
9 /* --------------------------------------------------------------------------------------- */
10 
11 int atmocor2(l2str *l2rec, aestr *aerec, int32_t ip) {
12  static int firstCall = 1;
13  static int want_nirRrs = 0;
14  static int want_nirLw = 0;
15  static int want_mumm = 0;
16  static int want_ramp = 1;
17  static int32_t aer_iter_max = 1;
18  static int32_t aer_iter_min = 1;
19 
20  static float pi = PI;
21  static float radeg = RADEG;
22  static float badval = BAD_FLT;
23  static float badchl = BAD_FLT;
24  static float p0 = STDPR;
25  static float df = 0.33;
26  static float cbot = 0.7;
27  static float ctop = 1.3;
28  static float seed_chl = 0.0;
29  static float seed_green = 0.0;
30  static float seed_red = 0.0;
31  static float nir_chg = 0.02;
32  static float glint_min = GLINT_MIN;
33  static float cslp;
34  static float cint;
35 
36  static int32_t green;
37  static int32_t red;
38  static int32_t nir_s;
39  static int32_t nir_l;
40  static int32_t swir_s;
41  static int32_t swir_l;
42  static int32_t aer_s;
43  static int32_t aer_l;
44  static int32_t daer;
45  static int32_t aer_base;
46  static int32_t nwvis;
47  static int32_t nwave;
48  static float *wave;
49  int32_t nWaveCovariance = 1;
50 
51  l1str *l1rec = l2rec->l1rec;
52  filehandle *l1file = l1rec->l1file;
53  uncertainty_t *uncertainty=l1rec->uncertainty;
54  int32_t proc_uncertainty = input->proc_uncertainty;
55 
56  int32_t sensorID = l1file->sensorID;
57  int32_t brdf_opt = input->brdf_opt;
58  int32_t aer_opt = input->aer_opt;
59  int32_t glint_opt = input->glint_opt;
60  int32_t cirrus_opt = input->cirrus_opt;
61 
62  float *Fo = l1rec->Fo;
63  float *Fobar = l1file->Fobar;
64  float fsol = l1rec->fsol;
65  float solz = l1rec->solz [ip];
66  float senz = l1rec->senz [ip];
67  float delphi = l1rec->delphi[ip];
68  float ws = l1rec->ws [ip];
69  float gc = l1rec->glint_coef[ip];
70 
71  int32_t *aermodmin = &l2rec->aermodmin[ip];
72  int32_t *aermodmax = &l2rec->aermodmax[ip];
73  float *aermodrat = &l2rec->aerratio [ip];
74  int32_t *aermodmin2 = &l2rec->aermodmin2[ip];
75  int32_t *aermodmax2 = &l2rec->aermodmax2[ip];
76  float *aermodrat2 = &l2rec->aerratio2 [ip];
77  float *eps = &l2rec->eps [ip];
78 
79  float *TLg = &l1rec->TLg [ip * l1file->nbands];
80  float *Lw = &l2rec->Lw [ip * l1file->nbands];
81  float *nLw = &l2rec->nLw [ip * l1file->nbands];
82  float *La = &l2rec->La [ip * l1file->nbands];
83  float *taua = &l2rec->taua [ip * l1file->nbands];
84  float *t_sol = &l1rec->t_sol[ip * l1file->nbands];
85  float *t_sen = &l1rec->t_sen[ip * l1file->nbands];
86  float *brdf = &l2rec->brdf [ip * l1file->nbands];
87  float *Rrs = &l2rec->Rrs [ip * l1file->nbands];
88 
89  float *Rrs_unc;
90 
91  float *tg_sol = &l1rec->tg_sol[ip * l1file->nbands];
92  float *tg= &l1rec->tg[ip * l1file->nbands];; // double way absorption
93 
94  float *dsensor=NULL, *dLr=NULL, *dtaua=NULL, *dtg_sol=NULL;
95  float *dtg_sen=NULL, *dt_sol=NULL, *dt_sen=NULL, *dvc=NULL, *last_dtaua=NULL;
96  float dLt;
97  int dvc_fail = 0;
98 
99  float last_dtaua_aer_l, *derv_taua_l, **derv_rhorc, *derv_rhow_l, *derv_rh;
100 
101  float *derv_Lg_taua=NULL; //derivative of TLg[wave] to corresponding taua[nir_l]
102  float *drhown_nir=NULL;
103  float *dchl;
104  float *covariance_matrix;
105  static float *corr_coef_rhot;
106  int dim;
107  float *F1, *F1_temp, *F2, **COV, **corr_nir;
108  static float *delta_Lt = NULL;
109 
110  float *taur,*radref;
111  float *tLw, *dtLw;
112  float *rhown_nir;
113  float *tLw_nir, *dtLw_nir, *last_tLw_nir, *dlast_tLw_nir;
114  int32_t ib, ipb,inir,iw, i, j, ix1, ix2;
115  int32_t status;
116  int32_t iter, iter_max, iter_min, last_iter, iter_reset;
117  float mu, mu0;
118  int32_t nneg;
119  float airmass;
120  float *Ltemp, *dLtemp;
121  float chl;
122  int want_glintcorr;
123  float refl_nir = 0, last_refl_nir = 0;
124  float tindx;
125  float Ka = 0.8; /* 1.375 um channel transmittance for water vapor (<=1) Gao et al. 1998 JGR */
126  float *mbac_w;
127  static int nbands_ac=2;
128  static int *acbands_index=NULL;
129 
130  if (firstCall == 1) {
131  firstCall = 0;
132  nwave = l1file->nbands;
133  if ((wave = (float *) calloc(nwave, sizeof (float))) == NULL) {
134  printf("-E- : Error allocating memory to wave\n");
135  exit(FATAL_ERROR);
136  }
137  for (ib = 0; ib < nwave; ib++) {
138  wave[ib] = l1file->fwave[ib];
139  }
140  if (sensorID == CZCS) {
141  want_ramp = 0;
142  seed_chl = 0.01;
143  seed_green = 0.003;
144  seed_red = 0.00125;
145  }
146  if(!input->acbands_index)
147  input->acbands_index=(int *)malloc(nbands_ac*sizeof(int));
148  acbands_index=input->acbands_index;
149  nwvis = rdsensorinfo(l1file->sensorID, l1_input->evalmask, "NbandsVIS", NULL);
150  nir_s = bindex_get(input->aer_wave_short);
151  nir_l = bindex_get(input->aer_wave_long);
152  aer_base=bindex_get(input->aer_wave_base);
153  if (nir_s < 0 || nir_l < 0) {
154  printf("Aerosol selection bands %d and %d not available for this sensor\n",
155  input->aer_wave_short, input->aer_wave_long);
156  exit(1);
157  }
158  if (nir_l < nir_s) {
159  printf("Invalid aerosol selection bands: long (%d) must be greater than short (%d).\n",
160  input->aer_wave_long, input->aer_wave_short);
161  exit(1);
162  }
163  if (wave[nir_s] < 600) {
164  printf("Aerosol selection band(s) must be greater than 600nm");
165  exit(1);
166  }
167 
168  aer_s = nir_s;
169  aer_l = nir_l;
170  daer = MAX(nir_l - nir_s, 1);
171  cslp = 1. / (ctop - cbot);
172  cint = -cslp * cbot;
173  printf("Aerosol selection bands %d and %d\n", l1file->iwave[aer_s], l1file->iwave[aer_l]);
174 
175  switch (aer_opt) {
176  case AERWANGNIR:
177  case AERRHNIR:
178  case FIXANGSTROMNIR:
179  case FIXMODPAIRNIR:
180  want_nirLw = 1;
181  aer_iter_min = 1;
182  aer_iter_max = input->aer_iter_max;
183  printf("NIR correction enabled.\n");
184  break;
185  case AERMUMM:
186  case AERRHMUMM:
187  want_mumm = 1;
188  aer_iter_min = 3;
189  aer_iter_max = input->aer_iter_max;
190  printf("MUMM correction enabled.\n");
191  break;
192  case AERWANGSWIR:
193  case AERRHSWIR:
194  want_nirLw = 1;
195  aer_iter_min = 1;
196  aer_iter_max = input->aer_iter_max;
197  swir_s = bindex_get(input->aer_swir_short);
198  swir_l = bindex_get(input->aer_swir_long);
199  if (swir_s < 0 || swir_l < 0) {
200  printf("Aerosol selection bands %d and %d not available for this sensor\n",
201  input->aer_swir_short, input->aer_swir_long);
202  exit(1);
203  }
204  printf("NIR/SWIR switching correction enabled.\n");
205  break;
206  case AERRHMSEPS:
207  want_nirLw = 1;
208  aer_iter_min = 1;
209  aer_iter_max = input->aer_iter_max;
210  input->nbands_ac=nbands_ac;
211  acbands_index[0]=windex(input->aer_wave_short,wave,nwave);
212  acbands_index[1]=windex(input->aer_wave_long,wave,nwave);
213  printf("NIR correction enabled --> for multi-scattering epsilon.\n");
214  break;
215  case AERRHSM:
216  want_nirLw = 1; //This needs to be turned on, but a new SWIR water model is needed for it to work
217  aer_iter_min = 1;
218  aer_iter_max = input->aer_iter_max;
219  daer=1;
220  nbands_ac=input->nbands_ac;
221 
222  if(input->mbac_wave[0]<input->aer_wave_short){
223  printf("%s line %d: the first mbac wavelength %d shouldn't be shorter than aer_wave_short %d \n", __FILE__, __LINE__,input->mbac_wave[0],input->aer_wave_short);
224  exit(1);
225  }
226  for(ib=0;ib<nbands_ac;ib++){
227  acbands_index[ib]=0;
228  for(iw=aer_s;iw<nwave;iw++){
229  if(input->mbac_wave[ib]==l1file->iwave[iw]){
230  acbands_index[ib]=1;
231  break;
232  }
233  }
234  if(!acbands_index[ib]){
235  printf("%s line %d: the band %d specified in mbac_wave doesn't exist \n", __FILE__, __LINE__,input->mbac_wave[ib]);
236  exit(1);
237  }
238  acbands_index[ib]=windex(input->mbac_wave[ib],wave,nwave);
239  }
240 
241  printf("NIR correction enabled --> for spectral matching.\n");
242  break;
243  default:
244  if (input->aer_rrs_short >= 0.0 && input->aer_rrs_long >= 0.0) {
245  want_nirRrs = 1;
246  aer_iter_min = 3;
247  aer_iter_max = input->aer_iter_max;
248  printf("NIR correction via input Rrs enabled.\n");
249  }
250  break;
251  }
252  if (input->aer_iter_max < 1)
253  want_nirLw = 0;
254  if ( want_nirLw || want_nirRrs) {
255  if ((red = bindex_get(670)) < 0) {
256  if ((red = bindex_get(680)) < 0)
257  if ((red = bindex_get(620)) < 0)
258  if ((red = bindex_get(765)) < 0)
259  if ((red = bindex_get(655)) < 0)
260  if ((red = bindex_get(664)) < 0) /* added for MSI */ {
261  printf("%s line %d: can't find red band\n", __FILE__, __LINE__);
262  exit(1);
263  }
264  }
265  if ((green = bindex_get(550)) < 0) {
266  if ((green = bindex_get(555)) < 0)
267  if ((green = bindex_get(560)) < 0)
268  if ((green = bindex_get(565)) < 0) {
269  printf("%s line %d: can't find green band\n", __FILE__, __LINE__);
270  exit(1);
271  }
272  }
273  }
274 
275  if(uncertainty){
276  if ((delta_Lt = (float *)calloc(nwave, sizeof(float))) == NULL) {
277  printf("-E- : Error allocating memory to delta_Lt\n");
278  exit(FATAL_ERROR);
279  }
280  corr_coef_rhot = uncertainty->corr_coef_rhot;
281  }
282  }
283 
284  if ((taur = (float *)calloc(nwave, sizeof(float))) == NULL) {
285  printf("-E- : Error allocating memory to taur\n");
286  exit(FATAL_ERROR);
287  }
288  if ((tLw = (float *) calloc(nwave, sizeof (float))) == NULL) {
289  printf("-E- : Error allocating memory to tLw\n");
290  exit(FATAL_ERROR);
291  }
292  if ((rhown_nir = (float *) calloc(nwave, sizeof (float))) == NULL) {
293  printf("-E- : Error allocating memory to rhown_nir\n");
294  exit(FATAL_ERROR);
295  }
296  if ((tLw_nir = (float *) calloc(nwave, sizeof (float))) == NULL) {
297  printf("-E- : Error allocating memory to tLw_nir\n");
298  exit(FATAL_ERROR);
299  }
300  if ((last_tLw_nir = (float *) calloc(nwave, sizeof (float))) == NULL) {
301  printf("-E- : Error allocating memory to last_tLw_nir\n");
302  exit(FATAL_ERROR);
303  }
304  if ((Ltemp = (float *) calloc(nwave, sizeof (float))) == NULL) {
305  printf("-E- : Error allocating memory to Ltemp\n");
306  exit(FATAL_ERROR);
307  }
308  if ((mbac_w = (float *) calloc(nwave, sizeof (float))) == NULL) {
309  printf("-E- : Error allocating memory to mbac_w\n");
310  exit(FATAL_ERROR);
311  }
312  if ((radref = (float *) calloc(nwave, sizeof (float))) == NULL) {
313  printf("-E- : Error allocating memory to radref\n");
314  exit(FATAL_ERROR);
315  }
316 
317  if (uncertainty) {
318  init_uncertainty(uncertainty,0);
319  derv_Lg_taua=uncertainty->derv_Lg_taua;
320  drhown_nir=uncertainty->drhown_nir;
321  dchl=&uncertainty->dchl;
322  Rrs_unc = &l2rec->Rrs_unc[ip * l1file->nbands];
323  uncertainty->dRrs=Rrs_unc;
324  if (proc_uncertainty == 2) {
325  covariance_matrix = &l2rec->covariance_matrix[ip * l1file->nbands * l1file->nbands];
326  uncertainty->covaraince_matrix = covariance_matrix;
327  } else
328  covariance_matrix = uncertainty->pixel_covariance;
329 
330  dsensor = &uncertainty->dsensor[ip * nwave];
331  dLr = &uncertainty->dLr[ip * nwave];
332  dtg_sol = &uncertainty->dtg_sol[ip * nwave];
333  dtg_sen = &uncertainty->dtg_sen[ip * nwave];
334  dt_sol = &uncertainty->dt_sol[ip * nwave];
335  dt_sen = &uncertainty->dt_sen[ip * nwave];
336  dvc = &uncertainty->dvc[ip * nwave];
337  //dbrdf = &uncertainty->dbrdf[ip * nwave];
338 
339  if ((last_dtaua = (float *) calloc(nwave, sizeof(float))) == NULL) {
340  printf("-E- : Error allocating memory to Lunc\n");
341  exit(FATAL_ERROR);
342  }
343  if ((dtaua = (float *) calloc(nwave, sizeof(float))) == NULL) {
344  printf("-E- : Error allocating memory to Lunc\n");
345  exit(FATAL_ERROR);
346  }
347 
348  if ((dLtemp = (float *) calloc(nwave, sizeof(float))) == NULL) {
349  printf("-E- : Error allocating memory to dLtemp\n");
350  exit(FATAL_ERROR);
351  }
352  if ((dlast_tLw_nir = (float *) calloc(nwave, sizeof(float))) == NULL) {
353 
354  printf("-E- : Error allocating memory to last_tLw_nir\n");
355  exit(FATAL_ERROR);
356  }
357  if ((dtLw = (float *) calloc(nwave, sizeof(float))) == NULL) {
358  printf("-E- : Error allocating memory to tLw\n");
359  exit(FATAL_ERROR);
360  }
361  if ((dtLw_nir = (float *) calloc(nwave, sizeof(float))) == NULL) {
362  printf("-E- : Error allocating memory to dtLw_nir\n");
363  exit(FATAL_ERROR);
364  }
365  if ((derv_taua_l = (float *)calloc(nwave, sizeof(float))) == NULL) {
366  printf("-E- : Error allocating memory to derv_taua_l\n");
367  exit(FATAL_ERROR);
368  }
369  if ((derv_rhow_l = (float *)calloc(nwave, sizeof(float))) == NULL) {
370  printf("-E- : Error allocating memory to derv_rhow_l\n");
371  exit(FATAL_ERROR);
372  }
373  if ((derv_rh = (float *)calloc(nwave, sizeof(float))) == NULL) {
374  printf("-E- : Error allocating memory to derv_rh\n");
375  exit(FATAL_ERROR);
376  }
377  if ((derv_rhorc = (float **)calloc(nwave, sizeof(float *))) == NULL) {
378  printf("-E- : Error allocating memory to derv_rhorc\n");
379  exit(FATAL_ERROR);
380  }
381  for (ib = 0; ib < nwave; ib++) {
382  if ((derv_rhorc[ib] = (float *)calloc(nbands_ac, sizeof(float))) == NULL) {
383  printf("-E- : Error allocating memory to derv_rhorc[%d]\n", ib);
384  exit(FATAL_ERROR);
385  }
386  }
387 
389 
390  dim = nbands_ac + 4;
391  F1 = (float *)malloc(dim * sizeof(float));
392  F1_temp = (float *)malloc(dim * sizeof(float));
393  F2 = (float *)malloc(dim * sizeof(float));
394  COV = (float **)malloc(dim * sizeof(float *));
395  for (ib = 0; ib < dim; ib++)
396  COV[ib] = (float *)malloc(dim * sizeof(float));
397 
398  for (ib = 0; ib < dim; ib++)
399  for (iw = 0; iw < dim; iw++)
400  COV[ib][iw] = 0.;
401 
402  corr_nir = (float **)malloc(nwave * sizeof(float *));
403  for (ib = 0; ib < nwave; ib++)
404  corr_nir[ib] = (float *)malloc(nbands_ac * sizeof(float));
405  for (ib = 0; ib < nwave; ib++)
406  for (iw = 0; iw < nbands_ac; iw++)
407  corr_nir[ib][iw] = 0;
408 
409  for (iw = 0; iw < nwave; iw++) {
410  for (ib = 0; ib < nbands_ac; ib++) {
411  if ((acbands_index[ib]) >= iw)
412  corr_nir[iw][ib] = corr_coef_rhot[iw * nwave + acbands_index[ib]];
413  else
414  corr_nir[iw][ib] = corr_coef_rhot[acbands_index[ib] * nwave + iw];
415  }
416  }
417  }
418 
419  /* Initialize output values. If any input radiances are negative, just return. */
420 
421  status = 0;
422  iter = 0;
423  nneg = 0;
424 
425  for (ib = 0; ib < nwave; ib++) {
426  ipb = ip * nwave + ib;
427 
428  //t_sol [ib] = 1.0; leave them as rayleigh only
429  //t_sen [ib] = 1.0;
430  La [ib] = badval;
431  tLw [ib] = badval;
432  Lw [ib] = badval;
433  nLw [ib] = badval;
434  taua [ib] = badval;
435  rhown_nir[ib]=0.;
436  if (glint_opt != 2) TLg [ib] = 0.0;
437  brdf [ib] = 1.0;
438  Rrs [ib] = badval;
439 
440  l2rec->eps[ip] = badval;
441 
442  if (l2rec->l1rec->Lt[ipb] <= 0.0)
443  if (wave[ib] < 1000) nneg++; /* don't test SWIR */
444 
445  if(aer_opt==AERRHSM)
446  mbac_w[ib]=0.0;
447  if (uncertainty) {
448  dtaua[ib] = 0.;
449  //derv_Lg_taua[ib] = 0.;
450  last_dtaua[ib] = 0;
451  // Rrs_unc[ib] = 0.;
452  dvc[ib] = dvc[ib] / tg[ib] / l1rec->polcor[ipb];
453  }
454  }
455 
456  /* If any expected channels are negative */
457  if (nneg > 0) {
458  free(taur);
459  free(tLw);
460  free(rhown_nir);
461  free(tLw_nir);
462  free(last_tLw_nir);
463  free(Ltemp);
464  free(mbac_w);
465  free(radref);
466 
467  if(uncertainty){
468  free(dtLw);
469  free(dtLw_nir);
470  free(dlast_tLw_nir);
471  free(dLtemp);
472  free(last_dtaua);
473  free(dtaua);
474  free(derv_taua_l);
475  free(derv_rhow_l);
476  free(derv_rh);
477  for (ib = 0; ib < nwave; ib++)
478  free(derv_rhorc[ib]);
479  free(derv_rhorc);
480 
481  free(F1);
482  free(F2);
483  free(F1_temp);
484  for (ib = 0; ib < dim; ib++)
485  free(COV[ib]);
486  free(COV);
487 
488  for (ib = 0; ib < nwave; ib++)
489  free(corr_nir[ib]);
490  free(corr_nir);
491  }
492 
493  status = 1;
494  return (status);
495  }
496 
497  mu0 = cos(solz / radeg);
498  mu = cos(senz / radeg);
499  airmass = 1.0 / mu0 + 1.0 / mu;
500 
501  /* Remove pre-computed atmospheric effects */
502  /* --------------------------------------- */
503  for (ib = 0; ib < nwave; ib++) {
504 
505  ipb = ip * nwave + ib;
506 
507  /* Pressure correct the Rayleigh optical thickness */
508  taur[ib] = l1rec->pr[ip] / p0 * l1file->Tau_r[ib];
509 
510  /* Copy TOA radiance to temp var, eliminate missing bands */
511  Ltemp[ib] = l2rec->l1rec->Lt[ipb];
512 
513  /* Correct for ozone absorption. We correct for inbound and outbound here, */
514  /* then we put the inbound back when computing Lw. */
515  Ltemp[ib] = Ltemp[ib] / tg[ib];
516 
517  /* Do Cirrus correction - subtract off cirrus reflectance from Ltemp */
518  /* Ka is the 1.375 um transmittance of water vapor above cirrus clouds */
519  /* Add cirrus_opt to input */
520  if (cirrus_opt) Ltemp[ib] -= l1rec->rho_cirrus[ip] / Ka * Fo[ib] * mu0 / pi;
521 
522  /* Apply polarization correction */
523  Ltemp[ib] /= l1rec->polcor[ipb];
524 
525  /* Remove whitecap radiance */
526  Ltemp[ib] -= l1rec->tLf[ipb];
527 
528  /* Subtract the Rayleigh contribution for this geometry. */
529  Ltemp[ib] -= l1rec->Lr[ipb];
530 
531  /* If selected, correct for Ding and Gordon O2 absorption for the aerosol component (replace O2 losses) */
532  if (input->oxaband_opt == 1) {
533  Ltemp[ib] /= l1rec->t_o2[ipb];
534  }
535 
536 
537  /* calculate error in Ltemp */
538  if (uncertainty) {
539  dLt = (1 - l2rec->l1rec->Lt[ipb] * uncertainty->derv_pol[ipb] / l1rec->polcor[ipb]);
540  dLt *= (1 / l1rec->tg_sol[ipb] / l1rec->tg_sen[ipb] / l1rec->polcor[ipb]);
541  dLt = pow(dLt * dsensor[ib], 2);
542 
543  dLt += pow( l2rec->l1rec->Lt[ipb] / l1rec->tg_sol[ipb] / pow(l1rec->tg_sen[ipb], 2) / l1rec->polcor[ipb] * dtg_sen[ib], 2);
544  dLt += pow( l2rec->l1rec->Lt[ipb] / l1rec->tg_sen[ipb] / pow(l1rec->tg_sol[ipb], 2) / l1rec->polcor[ipb] * dtg_sol[ib], 2);
545 
546  dLtemp[ib] = dLt + dLr[ib] * dLr[ib];
547  }
548  radref[ib]=pi/(Fo[ib] * mu0);
549  }
550 
551 
552  /* Compute BRDF correction, if not dependent on radiance */
553  if (brdf_opt < FOQMOREL && brdf_opt > NOBRDF)
554  ocbrdf(l2rec, ip, brdf_opt, wave, nwvis, solz, senz, delphi, ws, -1.0, NULL, NULL, brdf);
555 
556  /* Initialize iteration loop */
557  chl = seed_chl;
558  iter = 0;
559  last_iter = 0;
560  iter_max = aer_iter_max;
561  iter_min = aer_iter_min;
562  iter_reset = 0;
563  last_refl_nir = 100.;
564  want_glintcorr = 0;
565 
566  /* new glint_opt usage - a 2 will use the simple TLg from atmocor1 */
567  if (glint_opt == 1 && l1rec->glint_coef[ip] > glint_min) {
568  iter_max = MAX(2, iter_max);
569  iter_min = MAX(2, iter_min);
570  want_glintcorr = 1;
571  }
572  if (glint_opt == 2) {
573  want_glintcorr = 1;
574  }
575 
576  if (aer_opt == AERWANGSWIR || aer_opt == AERRHSWIR) {
577  tindx_shi(l2rec, ip, &tindx);
578  if (tindx >= 1.05) {
579  iter_max = 1;
580  aer_s = swir_s;
581  aer_l = swir_l;
582  want_nirLw = 0;
583  } else {
584  aer_s = nir_s;
585  aer_l = nir_l;
586  want_nirLw = 1;
587  }
588  daer = MAX(aer_l - aer_s, 1);
589  }
590 NIRSWIR:
591 
592  if (want_nirLw || want_nirRrs) {
593  for (ib = 0; ib < nwave; ib++) {
594  last_tLw_nir[ib] = 0.0;
595  tLw_nir[ib] = 0.0;
596  Rrs[ib] = 0.0;
597 
598  if(uncertainty){
599  dtLw_nir[ib] = 0.;
600  dlast_tLw_nir[ib] = 0.;
601  }
602  }
603  Rrs[green] = seed_green;
604  Rrs[red ] = seed_red;
605  }
606 
607 
608  /* -------------------------------------------------------------------- */
609  /* Begin iterations for aerosol with corrections for non-zero nLw(NIR) */
610  /* -------------------------------------------------------------------- */
611 
612  if (aer_opt==AERRHSM ) {//&& tLw_nir[aer_s]>0.
613  for(ib=0;ib<nbands_ac;ib++)
614  mbac_w[acbands_index[ib]]=1.;
615  }
616 
617  if (uncertainty) {
618  for (ib = 0; ib < nwave; ib++)
619  delta_Lt[ib] = sqrt(dLtemp[ib] + dvc[ib] * dvc[ib]);
620  }
621 
622  while (!last_iter) {
623  iter++;
624  status = 0;
625 
626  if (uncertainty)
627  last_dtaua_aer_l = dtaua[aer_l];
628 
629  /* Initialize tLw as surface + aerosol radiance */
630  for (ib = 0; ib < nwave; ib++) {
631  tLw[ib] = Ltemp[ib];
632  if (uncertainty){
633  dtLw[ib] = sqrt(dLtemp[ib]);
634  last_dtaua[ib] = dtaua[ib];
635  }
636  }
637 
638  /* Compute and subtract glint radiance */
639  if (want_glintcorr) {
640 
641  if (glint_opt == 1)
642  glint_rad(iter, nwave, aer_s, aer_l, gc, airmass, mu0, Fo, taur, taua, tLw, TLg, uncertainty);
643 
644  for (ib = 0; ib < nwave; ib++) {
645  tLw[ib] -= TLg[ib];
646  }
647  }
648 
649  /* Adjust for non-zero NIR water-leaving radiances using input Rrs */
650  if (want_nirRrs) {
651 
652  rhown_nir[aer_s] = pi * input->aer_rrs_short;
653  rhown_nir[aer_l] = pi * input->aer_rrs_long;
654 
655  for (ib = aer_s; ib <= aer_l; ib += daer) {
656 
657  /* Convert NIR reflectance to TOA W-L radiance */
658  tLw_nir[ib] = rhown_nir[ib] / pi * Fo[ib] * mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
659 
660  /* Avoid over-subtraction */
661  if (tLw_nir[ib] > tLw[ib] && tLw[ib] > 0.0)
662  tLw_nir[ib] = tLw[ib];
663 
664  /* Remove estimated NIR water-leaving radiance */
665  tLw[ib] -= tLw_nir[ib];
666  }
667  }
668 
669 
670  /* Adjust for non-zero NIR water-leaving radiances using MUMM */
671  if (want_mumm) {
672 
673  get_rhown_mumm(l2rec, ip, aer_s, aer_l, rhown_nir);
674 
675  for (ib = aer_s; ib <= aer_l; ib += daer) {
676 
677  /* Convert NIR reflectance to TOA W-L radiance */
678  tLw_nir[ib] = rhown_nir[ib] / pi * Fo[ib] * mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
679 
680  /* Remove estimated NIR water-leaving radiance */
681  tLw[ib] -= tLw_nir[ib];
682  }
683  }
684 
685 
686  /* Adjust for non-zero NIR water-leaving radiances using IOP model */
687  if (want_nirLw) {
688 
689  ipb = ip*nwave;
690  get_rhown_eval(input->fqfile, Rrs, wave, aer_s, aer_l, nwave, &l1rec->sw_a_avg[ipb], &l1rec->sw_bb_avg[ipb], chl, solz, senz, delphi, rhown_nir, l2rec,ip);
691 
692  for (ib = aer_s; ib <= aer_l; ib += daer) {
693 
694  /* Convert NIR reflectance to TOA W-L radiance */
695  tLw_nir[ib] = rhown_nir[ib] / pi * Fo[ib] * mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
696  if(uncertainty)
697  dtLw_nir[ib] = Fo[ib] * mu0 / (pi * brdf[ib]) *
698  sqrt(pow(rhown_nir[ib] * t_sol[ib] * dt_sen[ib], 2) +
699  pow(rhown_nir[ib] * t_sen[ib] * dt_sol[ib], 2) +
700  pow(t_sen[ib] * t_sol[ib] * drhown_nir[ib], 2));
701 
702  /* Iteration damping */
703  tLw_nir[ib] = (1.0 - df) * tLw_nir[ib] + df * last_tLw_nir[ib];
704 
705  if(uncertainty)
706  dtLw_nir[ib] = sqrt( pow((1. - df) * dtLw_nir[ib], 2) + pow(df * dlast_tLw_nir[ib], 2));
707 
708  /* Ramp-up ?*/
709  if (want_ramp) {
710  if (chl > 0.0 && chl <= cbot) {
711  tLw_nir[ib] = 0.0;
712  if(uncertainty)
713  dtLw_nir[ib] = 0.0;
714  }
715  else if ((chl > cbot) && (chl < ctop)) {
716  tLw_nir[ib] *= (cslp * chl + cint);
717  if(uncertainty)
718  dtLw_nir[ib] *= (cslp * chl + cint); //error in chl should be included in the next
719  }
720  }
721 
722  /* Remove estimated NIR water-leaving radiance */
723  tLw[ib] -= tLw_nir[ib];
724  }
725 
726  if(uncertainty && tLw_nir[aer_l] > 0.){
727  for (ib = 0; ib <nbands_ac; ib ++){
728  inir=acbands_index[ib];
729  uncertainty->ratio_rhow[ib] = tLw_nir[inir] / tLw_nir[aer_l] * (Fo[aer_l] / Fo[inir]);
730  }
731  }
732  }
733 
734  if (uncertainty) {
735  if (aer_opt == AERRHMSEPS && fabs(uncertainty->derv_Lg_taua[aer_s]) > 0.) {
736  uncertainty->derv_eps_taua_s = -1/ (Ltemp[aer_l] - TLg[aer_l] - tLw_nir[aer_l]) * derv_Lg_taua[aer_s];
737  uncertainty->derv_eps_taua_s *= (Fo[aer_l] / Fo[aer_s]);
738 
739  uncertainty->derv_eps_taua_l = (Ltemp[aer_s] - TLg[aer_s] - tLw_nir[aer_s]) /
740  pow(Ltemp[aer_l] - TLg[aer_l] - tLw_nir[aer_l], 2) *
741  derv_Lg_taua[aer_l];
742  uncertainty->derv_eps_taua_l *= (Fo[aer_l] / Fo[aer_s]);
743  uncertainty->derv_eps_taua_l+=uncertainty->derv_eps_taua_s;
744  }
745 
746  for (ib = 0; ib < nwave; ib++)
747  uncertainty->derv_Lg_taua[ib] *= radref[ib];
748  }
749 
750  /* Compute the aerosol contribution */
751  if (status == 0) {
752  if (aer_opt != AERNULL)
753  status = aerosol(l2rec, aer_opt, aerec, ip, wave, nwave, aer_s, aer_l, Fo, tLw,
754  La, t_sol, t_sen, eps, taua, aermodmin, aermodmax, aermodrat,
755  aermodmin2, aermodmax2, aermodrat2, mbac_w);
756  else {
757  for (ib = 0; ib < nwave; ib++) {
758 
759  ipb = ip * nwave + ib;
760 
761  t_sol [ib] = 1.0;
762  t_sen [ib] = 1.0;
763  La [ib] = 0.0;
764  taua [ib] = 0.0;
765  *eps = 1.0;
766  }
767  }
768  }
769 
770 
771  /* Subtract aerosols and normalize */
772  if (status == 0) {
773 
774  for (ib = 0; ib < nwave; ib++) {
775 
776  /* subtract aerosol and normalize */
777  tLw[ib] = tLw[ib] - La[ib];
778  Lw [ib] = tLw[ib] / t_sen[ib] * tg_sol[ib];
779  nLw[ib] = Lw [ib] / t_sol[ib] / tg_sol[ib] / mu0 / fsol * brdf[ib];
780 
781  if(uncertainty){
782 
783  dt_sen[ib] = pow(uncertainty->derv_tsen_taua_l[ib] * last_dtaua_aer_l, 2);
784 
785  for(inir=0;inir<nbands_ac;inir++){
786  i=acbands_index[inir];
787  dt_sen[ib] += pow(uncertainty->derv_tsen_rhorc[ib][inir] * sqrt(dLtemp[i]) * radref[i], 2);
788  }
789 
790  dt_sen[ib] += pow(uncertainty->derv_tsen_rhow_l[ib] * dtLw_nir[aer_l] * radref[aer_l], 2);
791  dt_sen[ib] += pow(uncertainty->derv_tsen_rh[ib] * uncertainty->drh[ip],2);
792 
793  /* vicarious calibration contribution */
794  for(inir=0;inir<nbands_ac;inir++){
795  i=acbands_index[inir];
796  dt_sen[ib] += pow(uncertainty->derv_tsen_rhorc[ib][inir] *radref[i]* dvc[i], 2);
797  }
798 
799  //TBD,only works for mseps right now, need to tune for mbac.
800  for(inir=0;inir<nbands_ac;inir++){
801  i=acbands_index[inir];
802  dt_sen[ib] +=2*corr_nir[i][nbands_ac - 1] *(uncertainty->derv_tsen_rhorc[ib][inir] * radref[i] * dvc[i] *uncertainty->derv_tsen_rhorc[ib][nbands_ac - 1] * radref[aer_l] * dvc[aer_l]);
803  }
804 
805  if(dt_sen[ib]<0)
806  dt_sen[ib]=0.;
807  dt_sen[ib] = sqrt(dt_sen[ib]);
808  /* */
809 
810  dt_sol[ib] = pow(uncertainty->derv_tsol_taua_l[ib] * last_dtaua_aer_l, 2);
811 
812  for(inir=0;inir<nbands_ac;inir++){
813  i=acbands_index[inir];
814  dt_sol[ib] += pow(uncertainty->derv_tsol_rhorc[ib][inir] * sqrt(dLtemp[i]) * radref[i], 2);
815  }
816 
817  dt_sol[ib] += pow(uncertainty->derv_tsol_rhow_l[ib] * dtLw_nir[aer_l] * radref[aer_l], 2);
818  dt_sol[ib] += pow(uncertainty->derv_tsol_rh[ib] * uncertainty->drh[ip], 2);
819 
820  /* vicarious calibration contribution */
821  for(inir=0;inir<nbands_ac;inir++){
822  i=acbands_index[inir];
823  dt_sol[ib] += pow(uncertainty->derv_tsol_rhorc[ib][inir] * radref[i]* dvc[i], 2);
824  }
825 
826  for(inir=0;inir<nbands_ac;inir++){
827  i=acbands_index[inir];
828  dt_sol[ib] +=
829  2 * corr_nir[i][nbands_ac - 1] *
830  (uncertainty->derv_tsol_rhorc[ib][inir] * radref[i] * dvc[i] *
831  uncertainty->derv_tsol_rhorc[ib][nbands_ac - 1] * radref[aer_l] * dvc[aer_l]);
832  }
833 
834  if(dt_sol[ib]<0)
835  dt_sol[ib]=0;
836  dt_sol[ib] = sqrt(dt_sol[ib]);
837  /* */
838 
839  derv_taua_l[ib] = -derv_Lg_taua[ib] / (t_sol[ib] * t_sen[ib]) * (1 / radref[ib]);
840  derv_taua_l[ib] += -uncertainty->derv_La_taua_l[ib] / (t_sol[ib] * t_sen[ib]);
841  derv_taua_l[ib] +=
842  (-tLw[ib] / (t_sol[ib] * t_sen[ib] * t_sen[ib]) * uncertainty->derv_tsen_taua_l[ib]);
843  derv_taua_l[ib] +=
844  (-tLw[ib] / (t_sol[ib] * t_sol[ib] * t_sen[ib]) * uncertainty->derv_tsol_taua_l[ib]);
845 
846  for (inir = 0; inir < nbands_ac; inir++) {
847  derv_rhorc[ib][inir] =
848  -uncertainty->derv_La_rhorc[ib][inir] / (t_sol[ib] * t_sen[ib]);
849  derv_rhorc[ib][inir] += (-tLw[ib] / (t_sol[ib] * t_sen[ib] * t_sen[ib]) *
850  uncertainty->derv_tsen_rhorc[ib][inir]);
851  derv_rhorc[ib][inir] += (-tLw[ib] / (t_sol[ib] * t_sol[ib] * t_sen[ib]) *
852  uncertainty->derv_tsol_rhorc[ib][inir]);
853  }
854 
855  derv_rhow_l[ib] = -uncertainty->derv_La_rhow_l[ib] / (t_sol[ib] * t_sen[ib]);
856  derv_rhow_l[ib] +=
857  (-tLw[ib] / (t_sol[ib] * t_sen[ib] * t_sen[ib]) * uncertainty->derv_tsen_rhow_l[ib]);
858  derv_rhow_l[ib] +=
859  (-tLw[ib] / (t_sol[ib] * t_sol[ib] * t_sen[ib]) * uncertainty->derv_tsol_rhow_l[ib]);
860 
861  derv_rh[ib] = -uncertainty->derv_La_rh[ib] / (t_sol[ib] * t_sen[ib]);
862  derv_rh[ib] +=
863  (-tLw[ib] / (t_sol[ib] * t_sen[ib] * t_sen[ib]) * uncertainty->derv_tsen_rh[ib]);
864  derv_rh[ib] +=
865  (-tLw[ib] / (t_sol[ib] * t_sol[ib] * t_sen[ib]) * uncertainty->derv_tsol_rh[ib]);
866 
867  /* end of vicarious calibration */
868 
869  /* end of using a different way to calculate dnLw */
870 
871  dtaua[ib] = pow(uncertainty->derv_taua_taua_l[ib] * last_dtaua_aer_l, 2);
872  for (inir = 0; inir < nbands_ac; inir++) {
873  i = acbands_index[inir];
874  dtaua[ib] +=
875  pow(uncertainty->derv_taua_rhorc[ib][inir] * sqrt(dLtemp[i]) * radref[i], 2);
876  }
877 
878  dtaua[ib] += pow(uncertainty->derv_taua_rhow_l[ib] * dtLw_nir[aer_l] * radref[aer_l], 2);
879  dtaua[ib] += pow(uncertainty->derv_taua_rh[ib] * uncertainty->drh[ip], 2);
880 
881  /* vicarious calibration contribution */
882  for (inir = 0; inir < nbands_ac; inir++) {
883  i = acbands_index[inir];
884  dtaua[ib] += pow(uncertainty->derv_taua_rhorc[ib][inir] * radref[i] * dvc[i], 2);
885  }
886  dtaua[ib] +=
887  2 * corr_nir[aer_s][nbands_ac - 1] *
888  (uncertainty->derv_taua_rhorc[ib][nbands_ac - 1] * radref[aer_l] * dvc[aer_l] *
889  uncertainty->derv_taua_rhorc[ib][0] * radref[aer_s] * dvc[aer_s]);
890 
891  /* */
892  if (dtaua[ib] < 0)
893  dtaua[ib] = 0;
894 
895  dtaua[ib] = sqrt(dtaua[ib]);
896  }
897  }
898 
899  /*calculate covariance matrix for Rrs */
900  /* cov(Rrs1, Rrs2)=F1.Cov(X1,X2).F2 */
901  /* F1: matrix of partial derivative of Rrs1 to X1(matrix) dimension: 1*dim */
902  /* F2: matrix of partial derivative of Rrs2 to X2(matrix) dimension: dim*1 */
903  /* COV(X1,X2):variance-covariance matrix of (X1,X2) dimension: dim*dim */
904  /* X1(rhorc1, rhorc_NIR, rh, rhow_l,taua_l) */
905  /* X2(rhorc2, rhorc_NIR, rh, rhow_l,taua_l) */
906 
907  if (uncertainty) {
908  for (ib = 0; ib < nwave; ib++) {
909  if (proc_uncertainty == 2)
910  nWaveCovariance = nwave - ib;
911  for (iw = ib; iw < ib + nWaveCovariance; iw++) {
912  /* calculate cov(Rrs[ib],Rrs[iw]) */
913 
914  F1[0] = 1 / t_sol[ib] / t_sen[ib] / radref[ib];
915  for (ix1 = 1; ix1 <= nbands_ac; ix1++)
916  F1[ix1] = derv_rhorc[ib][ix1 - 1];
917  F1[ix1] = derv_rh[ib];
918  F1[ix1 + 1] = derv_rhow_l[ib];
919  F1[ix1 + 2] = derv_taua_l[ib];
920 
921  F2[0] = 1 / t_sol[iw] / t_sen[iw] / radref[iw];
922  for (ix1 = 1; ix1 <= nbands_ac; ix1++)
923  F2[ix1] = derv_rhorc[iw][ix1 - 1];
924  F2[ix1] = derv_rh[iw];
925  F2[ix1 + 1] = derv_rhow_l[iw];
926  F2[ix1 + 2] = derv_taua_l[iw];
927 
928  COV[0][0] = corr_coef_rhot[ib * nwave + iw] * delta_Lt[ib] * radref[ib] *
929  delta_Lt[iw] * radref[iw];
930  for (i = 0; i < nbands_ac; i++) {
931  ix1 = acbands_index[i];
932  COV[0][i + 1] =
933  corr_nir[ib][i] * delta_Lt[ib] * radref[ib] * delta_Lt[ix1] * radref[ix1];
934  }
935  COV[0][nbands_ac + 1] = 0;
936  COV[0][nbands_ac + 2] = 0;
937  COV[0][nbands_ac + 3] = 0;
938 
939  for (i = 0; i < nbands_ac; i++) {
940  ix2 = acbands_index[i];
941 
942  COV[i + 1][0] =
943  corr_nir[iw][i] * delta_Lt[iw] * radref[iw] * delta_Lt[ix2] * radref[ix2];
944 
945  for (j = 0; j < nbands_ac; j++) {
946  ix1 = acbands_index[j];
947  COV[i + 1][j + 1] = corr_nir[ix2][j] * delta_Lt[ix2] * radref[ix2] *
948  delta_Lt[ix1] * radref[ix1];
949  }
950  COV[i + 1][nbands_ac + 1] = 0;
951  COV[i + 1][nbands_ac + 2] = 0;
952  COV[i + 1][nbands_ac + 3] = 0;
953  }
954 
955  for (ix2 = 1; ix2 < 4; ix2++)
956  for (ix1 = 0; ix1 < dim; ix1++)
957  COV[nbands_ac + ix2][ix1] = 0;
958 
959  COV[nbands_ac + 1][nbands_ac + 1] = pow(uncertainty->drh[ip], 2);
960  COV[nbands_ac + 2][nbands_ac + 2] = pow(dtLw_nir[aer_l] * radref[aer_l], 2);
961  COV[nbands_ac + 3][nbands_ac + 3] = pow(last_dtaua_aer_l, 2);
962 
963  for (ix1 = 0; ix1 < dim; ix1++) {
964  F1_temp[ix1] = 0.;
965  for (ix2 = 0; ix2 < dim; ix2++)
966  F1_temp[ix1] += F1[ix2] * COV[ix2][ix1];
967  }
968  covariance_matrix[ib * nwave + iw] = 0.;
969  for (ix1 = 0; ix1 < dim; ix1++)
970  covariance_matrix[ib * nwave + iw] += F1_temp[ix1] * F2[ix1];
971 
972  covariance_matrix[ib * nwave + iw] *=
973  (brdf[ib] * brdf[iw] / (mu0 * mu0) / (fsol * fsol));
974  }
975  }
976  }
977 
978  /* Compute new estimated chlorophyll */
979  if (want_nirLw) {
980  refl_nir = Rrs[red];
981  for (ib = aer_s; ib <= aer_l; ib += daer) {
982  last_tLw_nir[ib] = tLw_nir[ib];
983  if(uncertainty)
984  dlast_tLw_nir[ib] = dtLw_nir[ib];
985  }
986  for (ib = 0; ib < nwvis; ib++) {
987  Rrs[ib] = nLw[ib] / Fobar[ib];
988  if(uncertainty) {
989  for (iw = ib; iw < nwvis; iw++)
990  covariance_matrix[ib * nwave + iw] /= (Fobar[ib] * Fobar[iw]);
991  Rrs_unc[ib] = sqrt(covariance_matrix[ib * nwave + ib]);
992  }
993  }
994  chl = get_default_chl(l2rec, Rrs);
995 
996  // if we passed atmospheric correction but the spectral distribution of
997  // Rrs is bogus (chl failed), assume this is a turbid-water case and
998  // reseed iteration as if all 670 reflectance is from water.
999 
1000  if (chl == badchl && iter_reset == 0 && iter < iter_max) {
1001  chl = 10.0;
1002  Rrs[red] = 1.0 * (Ltemp[red] - TLg[red]) / t_sol[red] / tg_sol[red] / mu0 / fsol / Fobar[red];
1003  iter_reset = 1;
1004  if(uncertainty)
1005  *dchl=0.;
1006  }
1007 
1008  // if we already tried a reset, and still no convergence, force one last
1009  // pass with an assumption that all red radiance is water component, and
1010  // force iteration to end. this will be flagged as atmospheric correction
1011  // failure, but a qualitatively useful retrieval may still result.
1012 
1013  if (chl == badchl && iter_reset == 1 && iter < iter_max) {
1014  chl = 10.0;
1015  Rrs[red] = 1.0 * (Ltemp[red] - TLg[red]) / t_sol[red] / tg_sol[red] / mu0 / fsol / Fobar[red];
1016  iter = iter_max; // so iter will trigger maxiter flag and ATMFAIL
1017  iter_reset = 2;
1018  if(uncertainty)
1019  *dchl=0.;
1020  }
1021  }
1022 
1023  } else {
1024 
1025  /* Aerosol determination failure */
1026  for (ib = 0; ib < nwave; ib++) {
1027  La [ib] = badval;
1028  tLw[ib] = badval;
1029  Lw[ib] = badval;
1030  nLw[ib] = badval;
1031  Rrs[ib] = badval;
1032  }
1033  }
1034 
1035  /* If NIR/SWIR switching, do secondary test for turbidity and reset if needed */
1036  if (iter == 1 && (aer_opt == AERWANGSWIR || aer_opt == AERRHSWIR)) {
1037  if (tindx >= 1.05 && status == 0) {
1038  for (ib = 0; ib < nwvis; ib++) {
1039  Rrs[ib] = nLw[ib] / Fobar[ib];
1040  }
1041  chl = get_default_chl(l2rec, Rrs);
1042  //printf("Checking turbidity %d %f %f %f\n",ip,tindx,chl,nLw[nir_l]);
1043  if (chl > 0 && (chl <= 1.0 || nLw[nir_l] < 0.08)) {
1044  iter_max = aer_iter_max;
1045  aer_s = nir_s;
1046  aer_l = nir_l;
1047  daer = MAX(aer_l - aer_s, 1);
1048  want_nirLw = 1;
1049  //printf("Reverting to NIR %d %f %f %f\n",ip,tindx,chl,nLw[nir_l]);
1050  goto NIRSWIR;
1051  } else
1052  l1rec->flags[ip] |= TURBIDW;
1053  }
1054  }
1055 
1056 
1057  /* Shall we continue iterating */
1058  if (status != 0) {
1059  last_iter = 1;
1060  } else if (iter < iter_min) {
1061  last_iter = 0;
1062  } else if (want_nirLw && (fabs(refl_nir - last_refl_nir) < fabs(nir_chg * refl_nir) || refl_nir < 0.0)) {
1063  last_iter = 1;
1064  } else if (want_mumm || want_nirRrs) {
1065  last_iter = 1;
1066  } else if (iter > iter_max) {
1067  last_iter = 1;
1068  }
1069 
1070  last_refl_nir = refl_nir;
1071  if (aer_opt==AERRHSM) {//&& tLw_nir[aer_s]>0.
1072  for (ib = aer_s; ib <= aer_l; ib ++) {
1073  if(sensorID==MODISA && ib<=aer_base && iter>2 && mbac_w[ib]!=0.){
1074  mbac_w[ib] = (iter*1.0/iter_max);
1075  mbac_w[ib] = exp(-7*mbac_w[ib]);
1076  }
1077  }
1078  }
1079 
1080  } /* end of iteration loop */
1081 
1082  l2rec->num_iter[ip] = iter;
1083 
1084 
1085  /* If the atmospheric correction failed, we don't need to do more. */
1086  if (status != 0) {
1087  free(taur);
1088  free(tLw);
1089  free(rhown_nir);
1090  free(tLw_nir);
1091  free(last_tLw_nir);
1092  free(Ltemp);
1093  free(mbac_w);
1094  free(radref);
1095 
1096  if(uncertainty){
1097  free(dtLw);
1098  free(dtLw_nir);
1099  free(dlast_tLw_nir);
1100  free(dLtemp);
1101  free(last_dtaua);
1102  free(dtaua);
1103  free(derv_taua_l);
1104  free(derv_rhow_l);
1105  free(derv_rh);
1106  for (ib = 0; ib < nwave; ib++)
1107  free(derv_rhorc[ib]);
1108  free(derv_rhorc);
1109 
1110  free(F1);
1111  free(F2);
1112  free(F1_temp);
1113  for (ib = 0; ib < dim; ib++)
1114  free(COV[ib]);
1115  free(COV);
1116 
1117  for (ib = 0; ib < nwave; ib++)
1118  free(corr_nir[ib]);
1119  free(corr_nir);
1120  }
1121  return (status);
1122 
1123  }
1124 
1125  /* If we used a NIR Lw correction, we record the tLw as it was predicted. */
1126  if (want_nirLw || want_nirRrs || want_mumm) {
1127  for (ib = aer_s; ib <= aer_l; ib += daer) {
1128  tLw[ib] = tLw_nir[ib];
1129  Lw [ib] = tLw[ib] / t_sen[ib] * tg_sol[ib];
1130  nLw[ib] = Lw [ib] / t_sol[ib] / tg_sol[ib] / mu0 / fsol * brdf[ib];
1131  Rrs[ib] = nLw[ib] / Fobar[ib];
1132  }
1133  }
1134 
1135  /* Convert water-leaving radiances from band-averaged to band-centered. */
1136  /* Switch mean solar irradiance from band-averaged to band centered also.*/
1137  if (l1_input->outband_opt >= 2) {
1138  nlw_outband(l1_input->evalmask, sensorID, wave, nwave, Lw, nLw, &l2rec->outband_correction[ip*nwave]);
1139  Fobar = l1file->Fonom;
1140  }
1141 
1142  /* Compute f/Q correction and apply to nLw */
1143  if (brdf_opt >= FOQMOREL) {
1144  ocbrdf(l2rec, ip, brdf_opt, wave, nwvis, solz, senz, delphi, ws, -1.0, nLw, Fobar, brdf);
1145  for (ib = 0; ib < nwvis; ib++) {
1146  nLw[ib] = nLw[ib] * brdf[ib];
1147 
1148  if(uncertainty) {
1149  Rrs_unc[ib] *= brdf[ib];
1150  for (iw = ib; iw < nwave; iw++)
1151  covariance_matrix[ib * nwave + iw] *= (brdf[ib] * brdf[iw]);
1152  }
1153  }
1154  }
1155 
1156  /* Compute final Rrs */
1157  for (ib = 0; ib < nwave; ib++) {
1158  if (ib != aer_s && ib != aer_l) {
1159  Rrs[ib] = nLw[ib] / Fobar[ib];
1160  l2rec->Rrs[ip * nwave + ib] = Rrs[ib];
1161  }
1162  }
1163 
1164  if (uncertainty) {
1165  for (ib = nwvis; ib < nwave; ib++)
1166  for (iw = ib; iw < nwave; iw++)
1167  covariance_matrix[ib * nwave + iw] /= (Fobar[ib] * Fobar[iw]);
1168 
1169  for (ib = aer_s; ib < nwave; ib += daer) {
1170  if (Rrs_unc[ib] < 0)
1171  Rrs_unc[ib] = BAD_FLT;
1172  else {
1173  for (iw = ib; iw < nwave; iw++)
1174  covariance_matrix[ib * nwave + iw] *= (brdf[ib] * brdf[iw]);
1175  }
1176  }
1177  }
1178 
1179  /* Compute final chl from final nLw (needed for flagging) */
1180  l2rec->chl[ip] = get_default_chl(l2rec, Rrs);
1181 
1182  /* if there is no lower bounding for aerosol selection, the error can't be quantified */
1183  if (uncertainty) {
1184  if (*aermodmin == *aermodmax || *aermodmin2 == *aermodmax2 || dvc_fail) {
1185  for (ib = 0; ib < nwave; ib++) {
1186  Rrs_unc[ib] = BAD_FLT;
1187  Rrs[ib] = BAD_FLT;
1188  for (iw = 0; iw < nwave; iw++)
1189  covariance_matrix[ib * nwave + iw] = BAD_FLT;
1190  l2rec->chl[ip] = BAD_FLT;
1191  }
1192  } else {
1193  l2rec->chl_unc[ip] = *dchl;
1194  }
1195  }
1196 
1197  /*Determine Raman scattering contribution to Rrs*/
1198  run_raman_cor(l2rec, ip);
1199 
1200  free(taur);
1201  free(tLw);
1202  free(rhown_nir);
1203  free(tLw_nir);
1204  free(last_tLw_nir);
1205  free(Ltemp);
1206  free(mbac_w);
1207  free(radref);
1208 
1209  if(uncertainty){
1210  free(dtLw);
1211  free(dtLw_nir);
1212  free(dlast_tLw_nir);
1213  free(dLtemp);
1214  free(last_dtaua);
1215  free(dtaua);
1216  free(derv_taua_l);
1217  free(derv_rhow_l);
1218  free(derv_rh);
1219  for (ib = 0; ib < nwave; ib++)
1220  free(derv_rhorc[ib]);
1221  free(derv_rhorc);
1222 
1223  free(F1);
1224  free(F2);
1225  free(F1_temp);
1226  for (ib = 0; ib < dim; ib++)
1227  free(COV[ib]);
1228  free(COV);
1229 
1230  for (ib = 0; ib < nwave; ib++)
1231  free(corr_nir[ib]);
1232  free(corr_nir);
1233  }
1234 
1235  return (status);
1236 }
1237 
1238 
1239 
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
MAX
#define MAX(A, B)
Definition: swl0_utils.h:25
red
void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf, int ic1, int jc1, int jcf, int kc, float ***c, float **s)
j
int j
Definition: decode_rs.h:73
status
int status
Definition: l1_czcs_hdf.c:32
color_dtdb.gc
int gc
Definition: color_dtdb.py:186
F2
#define F2
Definition: anc.h:30
nlw_outband
void nlw_outband(int32_t evalmask, int32_t sensorID, float wave[], int32_t nwave, float Lw[], float nLw[], float outband_correction[])
Definition: nlw_outband.c:3
FIXMODPAIRNIR
#define FIXMODPAIRNIR
Definition: l12_parms.h:27
atrem_corl1.h
FOQMOREL
#define FOQMOREL
Definition: l12_parms.h:53
AERRHSM
#define AERRHSM
Definition: l12_parms.h:37
NULL
#define NULL
Definition: decode_rs.h:63
FIXANGSTROMNIR
#define FIXANGSTROMNIR
Definition: l12_parms.h:29
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
get_rhown_eval
void get_rhown_eval(char *fqfile, float Rrs[], float wave[], int32_t nir_s, int32_t nir_l, int32_t nwave, float aw[], float bbw[], float chl, float solz, float senz, float phi, float rhown[], l2str *l2rec, int32_t ip)
Definition: get_rhown_nir.c:286
tindx_shi
void tindx_shi(l2str *l2rec, int32_t ip, float *tindx)
Definition: turbid.c:4
atmocor2
int atmocor2(l2str *l2rec, aestr *aerec, int32_t ip)
Definition: atmocor2.c:11
pi
const double pi
Definition: VcstCmnConsts.h:473
mu
float mu
Definition: atrem_corl1v2.h:324
aerosol
int aerosol(l2str *l2rec, int32_t aer_opt_in, aestr *aerec, int32_t ip, float wave[], int32_t nwave, int32_t iwnir_s_in, int32_t iwnir_l_in, float F0_in[], float La1_in[], float La2_out[], float t_sol_out[], float t_sen_out[], float *eps, float taua_out[], int32_t *modmin, int32_t *modmax, float *modrat, int32_t *modmin2, int32_t *modmax2, float *modrat2, float *mbac_w)
Definition: aerosol.c:5651
eps
double eps
Definition: gha2000.c:3
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
AERRHNIR
#define AERRHNIR
Definition: l12_parms.h:24
hico.proc_hico.senz
senz
Definition: proc_hico.py:239
input
instr * input
Definition: main_l2binmatch.cpp:27
glint_rad
void glint_rad(int32_t iter, int32_t nband, int32_t nir_s, int32_t nir_l, float glint_coef, float airmass, float mu0, float F0[], float taur[], float taua[], float La[], float TLg[], uncertainty_t *uncertainty)
Definition: glint.c:40
l12_proto.h
PI
#define PI
Definition: l3_get_org.c:6
GLINT_MIN
#define GLINT_MIN
Definition: l1.h:60
TURBIDW
#define TURBIDW
Definition: l2_flags.h:22
bindex_get
int bindex_get(int32_t wave)
Definition: windex.c:45
NOBRDF
#define NOBRDF
Definition: l12_parms.h:50
tg_sol
float tg_sol[NBANDS]
Definition: atrem_corl1.h:127
get_rhown_mumm
void get_rhown_mumm(l2str *l2rec, int32_t ip, int32_t nir_s, int32_t nir_l, float rhown[])
Definition: mumm.c:74
AERRHSWIR
#define AERRHSWIR
Definition: l12_parms.h:33
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
create_images.Rrs
def Rrs
Definition: create_images.py:92
AERMUMM
#define AERMUMM
Definition: l12_parms.h:32
init_uncertainty
void init_uncertainty(uncertainty_t *uncertainty, int ifscan)
Definition: uncertainty.c:177
RADEG
#define RADEG
Definition: czcs_ctl_pt.c:5
AERRHMSEPS
#define AERRHMSEPS
Definition: l12_parms.h:36
STDPR
#define STDPR
Definition: l12_parms.h:85
run_raman_cor
void run_raman_cor(l2str *l2rec, int ip)
Definition: raman.c:817
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
get_default_chl
float get_default_chl(l2str *l2rec, float Rrs[])
Definition: get_chl.c:679
AERRHMUMM
#define AERRHMUMM
Definition: l12_parms.h:35
fabs
#define fabs(a)
Definition: misc.h:93
F1
#define F1
Definition: anc.h:29
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
CZCS
#define CZCS
Definition: sensorDefs.h:17
mu0
float mu0
Definition: atrem_corl1v2.h:324
rhown_nir
void rhown_nir(char *fqfile, float chl, float aw[], float bbw[], float Rrs[], float wave[], int32_t nwave, float solz, float senz, float phi, int32_t nir_s, int32_t nir_l, float rhown[], l2str *l2rec, int32_t ip)
Definition: get_rhown_nir.c:108
rdsensorinfo
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
AERWANGNIR
#define AERWANGNIR
Definition: l12_parms.h:25
i
int i
Definition: decode_rs.h:71
sensorID
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:91
ocbrdf
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
MODISA
#define MODISA
Definition: sensorDefs.h:19
AERWANGSWIR
#define AERWANGSWIR
Definition: l12_parms.h:31
AERNULL
#define AERNULL
Definition: l12_parms.h:39


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